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Seven  models  for  computing  underwater  radiances  and  irradiances  by  numerical  solution  of  the  radiative 
transfer  equation  are  compared.  The  models  are  applied  to  the  solution  of  several  problems  drawn  from 
optical  oceanography.  The  problems  include  highly  absorbing  and  highly  scattering  waters,  scattering 
by  molecules  and  by  particulates,  stratified  water,  atmospheric  effects,  surface-wave  effects,  bottom 
effects,  and  Raman  scattering.  The  models  provide  consistent  output,  with  errors  (resulting  from  Monte 
Carlo  statistical  fluctuations)  in  computed  irradiances  that  are  seldom  larger,  and  are  usually  smaller, 
than  the  experimental  errors  made  in  measuring  irradiances  when  using  current  oceanographic 
instrumentation.  Computed  radiances  display  somewhat  larger  errors.  > 


1 .  Introduction 

Various  numerical  models  are  now  in  use  for  comput¬ 
ing  underwater  irradiances  and  radiance  distributions. 
These  models  were  designed  to  address  a  wide  range 
of  oceanographic  problems.  The  models  are  based 
on  various  simplifying  assumptions,  have  differing 
levels  of  sophistication  in  their  representation  of 
physical  processes,  and  use  several  different  numeri¬ 
cal  solution  techniques. 

In  spite  of  the  increasingly  important  roles  these 
numerical  models  are  playing  in  optical  oceanogra¬ 
phy,  the  models  remain  incompletely  validated  in  the 
sense  that  their  outputs  have  not  been  extensively 
compared  with  measured  values  of  the  quantities 
they  predict.  This  desirable  model-data  comparison 
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is  not  presently  possible  because  the  requisite  compre¬ 
hensive  oceanic  optical  data  sets  are  not  available. 
Such  data  sets  must  contain  simultaneous  measure¬ 
ments  of  the  inherent  optical  properties  of  the  sea 
water  (e.g.,  the  absorption  and  scattering  coefficients 
and  the  scattering  phase  function),  environmental 
parameters  (e.g.,  the  sky  radiance  distribution  and 
sea  state),  and  radiometric  quantities  (e.g.,  the  com¬ 
plete  radiance  distribution  or  various  irradiances). 

The  inherent  optical  properties  and  the  environmen¬ 
tal  parameters  are  needed  as  input  to  the  numerical 
models;  the  radiometric  variables  are  the  quantities 
predicted  by  the  models.  Current  developments  in 
oceanic  optical  instrumentation  and  measurement 
methodologies  give  cause  for  hope  that  data  sets  that 
are  adequate  for  comprehensive  model-data  compari¬ 
sons  will  become  available  within  the  next  few  years. 

Meanwhile,  our  faith  in  these  models’  predictions 
rests  on  careful  debugging  of  computer  codes,  inter¬ 
nal  checks  such  as  conservation  of  energy  or  known 
relations  between  inherent  and  apparent  optical  prop¬ 
erties,  simulation  of  a  few  grossly  simplified  situa¬ 
tions  for  which  analytical  solutions  of  the  radiative 
transfer  equation  are  available,  and  comparison  (some¬ 
times  indirect)  with  incomplete  data  sets.  An  addi¬ 
tional  worthwhile  check  on  the  various  models  can  be 
made  by  applying  them  to  a  common  set  of  realistic 
problems.  Such  model-model  comparisons  help  to 
identify  errors  in  coding  or  weaknesses  in  the  math¬ 
ematical  representation  of  physical  phenomena,  quan¬ 
tify  numerical  errors  particular  to  the  various  solu¬ 
tion  algorithms,  determine  optimum  numerical 
techniques  for  simulation  of  particular  physical  phe- 
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nomena,  and  determine  which  models  might  be  appro¬ 
priate  for  inclusion  in  a  future  library  of  underwater 
radiative  transfer  codes  corresponding  to  those  now 
available  for  atmospheric  radiative  transfer  modeling 
(such  as  LOWTRAN1). 

In  March  1991  the  Oceanic  Optics  Program  of  the 
Office  of  Naval  Research  sponsored  a  workshop  to 
foster  a  close  examination  of  the  various  models  now 
in  use  and  to  begin  the  process  of  model-model 
comparison.  This  paper  reports  the  results  of  that 
comparison.  The  models  being  evaluated  are  de¬ 
scribed  in  Section  2.  During  the  workshop  the  par¬ 
ticipants  defined  a  set  of  canonical  (standard)  prob¬ 
lems  for  use  in  model  comparisons.  These  problems 
are  documented  in  Section  3.  Section  4  presents 
selected  results  obtained  when  the  models  of  Section 
2  are  applied  to  the  problems  of  Section  3. 

2.  Numerical  Models 

All  the  numerical  models  compared  here  generate  an 
approximate  solution  to  the  time-independent,  mono¬ 
chromatic  radiative  transfer  equation  in  one  spatial 
dimension: 

dL(v,  m  <|>)  r<  i%  ,  ,  ff  ,,  , 

M- - - =  -L(v,  H,  4>)  +  w0(t)  I  I  L( t;  M-  ,  <b  ) 

X  P(t;  pL'»  <K  -*■  4>)dpL'd<f»' 

+  S(t;  jt,  <}>).  (1) 

Here  L(t;  p.,  <J>)  is  the  unpolarized  spectral  radiance 
(at  wavelength  A,  omitted  for  brevity)  at  optical  depth 
t  and  in  direction  (p.,  <J>),  «*>o  is  the  scattering-to- 
attenuation  ratio,  (3  is  the  scattering  phase  function, 
and  S  represents  any  internal  source  of  radiance. 
The  depth  t  is  measured  positive  downward  from  the 
mean  sea  surface,  and  the  polar  angle  8  =  cos-1  p.  is 
measured  from  the  nadir  direction.  (See  Table  1  for 
a  list  of  symbols  and  their  units  and  definitions.) 
To  solve  Eq.  (1)  within  a  water  body,  it  is  necessary  to 
specify  (a)  the  inherent  optical  properties  of  the  water 
body,  <D0  and  (J;  (b)  the  distribution  of  internal  sources 
S;  (c)  the  radiance  distribution  that  is  externally 
incident  upon  the  boundaries  of  the  water  body;  and 
(d)  the  physical  nature  of  the  boundaries  themselves. 

The  models  differ  primarily  in  the  mathematical 
techniques  used  to  solve  Eq.  (1)  and  in  the  treatment 
of  boundary  conditions  at  the  sea  surface.  Two  of 
the  models  described  below  (models  II  and  DO)  use 
analytical  (invariant  imbedding  and  discrete  ordi¬ 
nates)  techniques  for  solving  Eq.  (1),  and  five  of  the 
models  (MC1-MC5)  use  probabilistic  (Monte  Carlo) 
techniques.  Each  of  the  models,  as  applied  to  the 
solution  of  the  canonical  problems  defined  in  Section 

3,  solves  Eq.  (1)  for  a  plane-parallel  water  body  that  is 
laterally  homogeneous  but  may  be  inhomogeneous 
with  depth.  The  upper  boundary  of  the  water  body 
is  the  windblown,  random  air-water  interface.  The 
lower  boundary  is  either  an  infinitely  thick  layer  of 
water  below  the  greatest  depth  of  interest  or  an 


Tablal.  Significant  Symbols,  Unit*,  and  Definition* 


Symbol 

Units 

Definition 

run 

Wavelength  of  tight 

* 

deg 

Scattering  angle;  0  s  iti  180° 

e 

deg 

Polar  angle  of  photon  travel,  measured 
from  the  nadir,  0  <,  8  <  180° 

<t> 

deg 

Azimuthal  angle  of  photon  travel,  mea¬ 
sured  counterclockwise  (looking  down¬ 
ward)  from  the  downwind  direction, 

0  s  d>  <  360° 

— 

p.  s  cos  8;  alternate  way  to  specify  the 
polar  angle,  -1  <  p.  <  1 

9* 

deg 

Polar  viewing  angle;  8t.  =  180°  -  0 

<t>„ 

deg 

Azimuthal  viewing  angle;  <)>„  =  180°  +  <t> 

ft 

sr 

Solid  angle;  a  differential  element  is  dft  = 
sin  8d0d<b  =  dpd<|> 

— 

Set  of  all  downward  directions; } a  dft  = 
2irsr 

w 

— u 

— 

Set  of  all  upward  directions;  j-  dft  =  2ir 

Z 

m 

Geometric  depth,  measured  positive 
downward 

T 

— 

Optical  depth,  measured  positive  down¬ 
ward:  T  =  I*  c(z)dz 

ff 

— 

Standard  deviation  of  the  surface- wave 
slope 

Clw 

m-1 

Absorption  coefficient  for  pure  sea  water 

aP 

m"1 

Absorption  coefficient  for  suspended  par¬ 
ticles 

a 

m“l 

Total  absorption  coefficient:  a  =  aw  +  ap 

bw 

m~l 

Scattering  coefficient  for  pure  sea  water 

bP 

m~l 

Scattering  coefficient  for  suspended  par¬ 
ticles 

b 

m_1 

Total  scattering  coefficient:  b  =  bu,  +  bp 

c 

m_l 

Total  attenuation  coefficient:  c  s  a  +  b 

p 

nr1  sr"1 

Volume  scattering  function 

p 

sr-1 

Scattering  phase  function,  p  s  p/6 

u>0 

— 

Scattering- to-attenuation  ratio,  o>o  =  b/c 

L 

Wm"2®"1 

Radiance  distribution,  L  =  L(z,  0,  <t>)  or 

nm"1 

£(t,  p,  4>) 

Lu 

W  m'2sr*1 

Nadir  radiance,  Lu  m  LJ t)  e  L(t,  0  =  180, 

nm'1 

4>  =  0) 

lm 

W  m'2sr"1 
nm-1 

Asymptotic  radiance  distribution 

kfj) 

— 

Asymptotic  diffuse  attenuation  coefficient 

s 

W  m~2sr_1 
nm-1 

Internal  source  of  radiance 

W  m"2  nm-1 

Irradiance  on  a  surface  perpendicular  to 
the  sun’s  rays 

Ed 

W  m~2  nm"1 

Downwelling  plane  irradiance:  Ed{-r)  = 
it,  4>)||t|dft 

Eu 

W  m"2  nm"1 

Upwelling plane  irradiance:  £„( t)  = 

J=  £(t,  p,  <f>)  1  it  |  dft 

Eou 

W  m"2  nm"1 

Upwelling  scalar  irradiance:  Em(j)  s 
|g_  L(t,  |t,  4>)dft 

opaque,  reflecting  bottom  at  a  finite  depth.  The 
models  all  assume  that  externally  applied  radiance  is 
incident  downward  upon  the  upper  side  of  the  air- 
water  surface.  The  models  are  all  monochromatic, 
and  there  are  no  internal  sources  of  radiance  (such  as 
bioluminescence).  However,  some  of  the  models  can 
simulate  inelastic-scattering  processes  by  sequential 
solutions  of  Eq.  (1).  For  example,  the  model  is  first 


20  December  1993  /  Vol.  32,  No.  36  /  APPLIED  OPTICS  7485 


¥ 


run  at  the  wavelength  of  excitation,  Aex,  to  compute 
the  energy  shifted  by  inelastic  scattering  from  Aex  to 
another  wavelength  A,  and  then  the  model  is  run 
again  at  A,  with  the  radiance  shifted  from  Aex  appear¬ 
ing  as  a  source  term  S  at  A.  A  particular  example  of 
S  used  in  this  treatment  of  Raman  scattering  is  given 
in  Appendix  A.  The  models  all  account  for  multiple 
scattering  and  can  use  realistic  scattering  phase 
functions  that  are  highly  peaked  in  forward  direc¬ 
tions,  as  is  the  case  for  sea  water. 

Several  of  the  models  have  additional  capabilities, 
such  as  the  computation  of  polarized  radiance  in  the 
Stokes  vector  format  and  the  simulation  of  azimuth- 
ally  anisotropic  random  air-water  surfaces.  These 
capabilities  are  not  evaluated  in  this  paper. 

All  but  one  of  the  models  directionally  discretizes 
Eq.  (1)  by  partitioning  the  set  of  all  directions,  5,  into 
a  grid  of  quadrilateral  regions  bounded  by  lines  of 
constant  p.  and  constant  <t>,  plus  two  polar  caps 
(collectively  called  quads).  The  fundamental  quan¬ 
tity  computed  by  these  models  is  the  quad-averaged 
radiance  defined  by 

L( t;  u,  y)  s  — -  I  j  L(t;  p,  d>)dpdd>.  (2) 

L(r,  u,  v)  is  physically  interpreted  as  the  average 
radiance  over  the  set  of  directions  (p,  4>)  contained  in 
the  uoth  quad,  Quv  (u  labels  p  bands  and  v  labels  4> 
bands),  which  subtends  a  solid  angle  of  size  In 

the  model  comparison  we  chose  to  use  24  <J>  bands  of 
uniform  width  A<f>  =  15°  and  20  p  bands  of  size  Ap  = 
0.1.  However,  a  polar  cap  with  Ap  =  0.1  has  a 
half-angle  of  0  =  25.8°,  which  is  much  larger  than  one 
would  normally  use  in  computing  the  nadir  or  zenith 
radiances.  Therefore  some  models  were  run  with  a 
slightly  different  p  spacing  and  smaller  polar  caps. 
The  remaining  model  (DO)  computes  the  radiance 
L(t;  p,  4>)  in  particular  (p,  <J>)  directions. 

We  now  briefly  describe  the  distinguishing  features 
of  the  various  models. 

A.  Model  II  [Invariant  Imbedding,  (Mobley)] 

The  integral  operator  of  Eq.  (2),  which  averages  any 
quantity  over  the  set  of  directions  (p,  <J>)  E  Quv,  can  be 
applied  to  Eq.  (1).  The  result  is  a  quad-averaged 
radiative  transfer  equation  in  which  L( t;  p,  4>)  is 
replaced  by  L( t;  u,v),  integration  over  all  directions  is 
replaced  by  summation  over  all  quads,  and  the  phase 
function  P(t;  p’,  4>'_— ►  p,  4>)  is  replaced  by  a  quad- 
averaged  quantity  0(t;  r,  s  -»  u,  v)  that  specifies  how 
much  of  the  radiance  initially  headed  into  quad  Qn 
gets  scattered  into  quad  Quv.  By  standard  tech¬ 
niques  of  Fourier  analysis  and  invariant  imbedding 
theory,  the  equations  for  the  L(t;  u,  v)  are  trans¬ 
formed  into  a  set  of  Riccati  differential  equations 
governing  the  depth  dependence  of  certain  reflec¬ 
tance  and  transmittance  functions  within  the  water 
body.  Depth  integration  of  the  Riccati  equations  (by 
a  high-order  Runge-Kutta  algorithm)  and  incorpora¬ 
tion  of  the  boundary  conditions  at  the  sea  surface  and 


bottom  lead  eventually  to  the  desired  L(r,  u,  v)  at  all 
depths.  These  mathematical  operations  are  out¬ 
lined  in  Mobley2  and  are  described  in  full  by  Mobley 
and  Preisendorfer.3  The  inherent  optical  properties 
of  the  water  body  can  vary  arbitrarily  with  depth. 
Absorption  and  scattering  are  built  up  as  sums  of 
terms  representing  the  contributions  by  pure  water, 
particles  of  various  types,  and  dissolved  substances. 

This  model  uses  a  Monte  Carlo  simulation  of  the 
windblown  sea  surface  to  evaluate  certain  quad- 
averaged,  bidirectional  reflectance  and  transmittance 
functions  that  describe  how  the  sea  surface  reflects 
and  transmits  radiance  that  is  incident  upon  the 
surface  from  above  and  below.  In  this  simulation, 
the  sea  surface  is  resolved  into  a  grid  of  triangular 
wave  facets  with  vertex  elevations  that  are  randomly 
determined  from  any  chosen  wave  slope-wind  speed 
spectrum  in  a  manner  described  by  Mobley  and 
Preisendorfer3  and  by  Preisendorfer  and  Mobley.4 
The  surface  simulation  allows  for  multiple  reflections 
of  rays  by  wave  facets  and  for  the  possibility  of 
shadowing  of  one  facet  by  another.  The  probabilis¬ 
tic  ray-tracing  calculations  for  setting  up  the  surface 
boundary  conditions  are  independent  of  the  analyti¬ 
cal  computations  within  the  water  body.  Moreover, 
because  the  ray  tracing  involves  only  the  surface- 
wave  facets,  for  which  it  is  assumed  that  there  is  no 
absorption,  no  rays  are  lost  to  absorption.  It  is 
therefore  computationally  feasible  to  trace  a  suffi¬ 
cient  number  of  rays  to  reduce  the  Monte  Carlo 
fluctuations  in  the  computed  bidirectional  surface 
functions  to  a  negligible  level. 

This  model  does  not  include  an  atmosphere  per  se. 
The  sky  radiance  that  is  incident  upon  the  sea  surface 
is  obtained  either  from  an  analytic  model  (e.g.,  a 
cardioidal  distribution  or  the  empirical  model  of 
Harrison  and  Coombes6)  or  from  the  output  of  a 
separately  run  atmospheric  radiative  transfer  model. 
In  the  simulation  of  problem  4,  below,  lowtran-7  was 
run  to  generate  the  sky  radiance  at  the  center  of  each 
of  the  (x — d>  quads;  that  value  was  then  taken  as  the 
average  sky  radiance  over  the  quad. 

The  bottom  boundary  can  be  either  an  infinitely 
thick  homogeneous  layer  of  water  below  some  depth 
Tmax  or  an  opaque  bottom  at  Tmax.  In  the  infinite- 
depth  case,  the  bidirectional  radiance  reflectance 
properties  of  the  infinite  layer  below  Tmax  are  obtained 
from  an  eigenmatrix  analysis  described  by  Preisendor¬ 
fer.6  The  same  analysis  yields  the  asymptotic  diffuse 
attenuation  coefficient,  and  the  asymptotic  radi¬ 
ance  distribution,  L*(p),  that  are  appropriate  for  the 
homogeneous  layer.  In  the  opaque-bottom  case,  the 
reflectance  properties  of  the  bottom  are  explicitly 
specified,  for  example,  as  a  Lambertian  surface  with  a 
given  irradiance  reflectance. 

The  chief  advantage  of  this  model  is  computational 
efficiency.  Solution  of  the  Riccati  differential  equa¬ 
tions  for  L  is  an  analytic  process,  and  thus  there  are 
no  Monte  Carlo  fluctuations  in  the  computed  radi¬ 
ances  (except  for  a  negligible  amount  introduced  by 
the  simulation  of  the  sea  surface).  In  particular, 
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both  upwelling  and  downwelling  radiances  are  com¬ 
puted  with  the  same  accuracy.  Moreover,  computa¬ 
tion  time  is  a  linear  function  of  depth,  so  that 
accurate  radiance  distributions  are  easily  obtained  at 
great  depths  (t  >  10).  Computation  time  depends 
only  mildly  on  quantities  such  as  the  scattering-to- 
attenuation  ratio,  surface  boundary  conditions,  and 
water  stratification.  The  associated  computer  code 
is  available  and  is  documented  by  Mobley.7 

8.  Model  DO  [Discrete  Ordinates  (Jin  and  Stamnes)] 

This  model  solves  Eq.  (1)  directly  without  applying 
the  quad-averaging  defined  by  Eq.  (2).  The  radiance 
is  expanded  into  a  Fourier  cosine  series,  L(r,  g,  4>)  = 
Lm( t,  p)cos(4>  -  4>o),  and  the  phase  function  is 
expanded  into  a  series  of  2 N  Legendre  polynomials, 

P(t;  p',  4>'  —  p,  4>)  s  0(r;  cos  i|») 

21V- 1 

=  2  (2/  +  1)£/(t)P/(cos  v|»), 

1=0 

where  gi(t)  is  the  expansion  coefficient  and  is  the 
scattering  angle.  The  advantage  of  these  expansions 
is  that  the  azimuthal  dependence  is  isolated,  in  the 
sense  that  2 N  independent  equations  for  the  Fourier 
coefficients  Lm( t,  p)  are  obtained: 

Hr.m( t  p)  C1 

P - ^ - =  -Lm(r,  p)  +  0»o(t)  J  Lm( T,  p') 

x  pm(T;  p',  p)dp'  +  Sm{ T,  p), 

where 

j  2N-1 

Pm(r;  p',  p)  =  o  2  +  ijgM 

Z  l=m 

{l  -  m)l 

X  WTmfF'"MPrM- 

Here  P^p)  is  the  associated  Legendre  polynomial. 

The  atmosphere  and  the  ocean  are  divided  into  a 
suitable  number  of  layers  to  adequately  resolve  the 
optical  properties  of  each  of  the  two  media.  Each 
layer  is  taken  to  be  homogeneous,  but  the  optical 
properties  are  allowed  to  vary  from  layer  to  layer. 
For  a  homogeneous  medium,  only  one  layer  is  re¬ 
quired.  At  the  interface  between  the  ocean  and  the 
atmosphere  (assumed  to  be  flat),  Fresnel’s  formula  is 
used  to  compute  the  appropriate  reflection  and  trans¬ 
mission  coefficients,  and  Snell’s  law  is  applied  to 
account  for  the  refraction  taking  place  there. 

The  integral  term  in  each  of  these  azimuth- 
independent  equations  is  then  approximated  by  a 
Gaussian  quadrature  sum  with  2 A/j  terms  (streams) 
in  the  atmosphere  and  2 N2  terms  in  the  ocean,  so  that 
there  are  21Vj  streams  in  the  refractive  region  of 
ocean  that  communicate  directly  with  the  atmo¬ 
sphere  and  2N2  -  2N)  streams  in  the  total  reflection 
region  of  the  ocean.  In  this  way  the  integro- 
differential  equation  is  transformed  into  a  system  of 


coupled  ordinary  differential  equations  that  is  solved 
by  the  discrete  ordinate  method,  as  described  in  more 
detail  elsewhere,8  subject  to  appropriate  boundary 
conditions  at  the  top  of  the  atmosphere  and  the 
bottom  of  the  ocean.  The  basic  discrete-ordinate 
method  used  here  is  described  and  thoroughly  docu¬ 
mented  in  previous  publications.9-11  The  modifica¬ 
tions  required  to  apply  the  method  to  a  system 
consisting  of  two  adjacent  media  with  different  indi¬ 
ces  of  refraction  are  described  by  Jin  and  Stamnes.8 

This  method  has  the  following  unique  features: 
(i)  Because  the  solution  is  analytic,  the  computational 
speed  is  completely  independent  of  individual  layer 
and  total  optical  thickness,  which  may  be  taken  to  be 
arbitrarily  large.  The  computational  speed  is  di¬ 
rectly  proportional  to  the  number  of  horizontal  layers 
used  to  resolve  the  optical  properties  in  the  atmo¬ 
sphere  and  ocean,  (ii)  Accurate  irradiances  are  ob¬ 
tained  with  just  a  few  streams,  which  makes  the  code 
very  efficient,  (iii)  Because  the  solution  is  analytic, 
radiances  and  irradiances  can  be  returned  at  arbi¬ 
trary  optical  depths  unrelated  to  the  computational 
levels,  (iv)  The  DO  method  is  essentially  a  matrix 
eigenvalue-eigenvector  solution,  from  which  the  as¬ 
ymptotic  solution  is  automatically  obtained.  The 
smallest  eigenvalue  is  k *,  and  the  associated  eigenvec¬ 
tor  is  Lx. 

Desirable  and  possible  extensions  of  the  method 
include  (i)  the  computation  of  inelastic-scattering 
effects  to  treat  phenomena  such  as  Raman  scattering 
and  (ii)  the  inclusion  of  a  windblown  surface  to 
simulate  the  basic  features  of  sea-surface  roughness. 
These  extensions  would  require  some  modifications 
of  the  existing  computer  code. 

C.  Model  MCI  [Monte  Carlo  (Gordon)J 

This  model  simulates  radiative  transfer  in  both  the 
ocean  and  the  atmosphere,  as  coupled  across  a  wind- 
roughened  interface.  The  code  is  designed  to  simu¬ 
late  irradiances  as  a  function  of  depth  for  computa¬ 
tion  of  the  irradiance  reflectance,  Eu/Ed,  and  diffuse 
attenuation  functions  such  as  Kd  =  ~&{\nEd)/dz. 
The  nadir- viewing  radiance,  Lu,  is  also  computed  as  a 
function  of  depth  for  the  computation  of  Q  =  Eu/Lu. 
The  optical  properties  of  the  ocean  are  continuously 
stratified  in  the  vertical.  They  can  be  specified  as 
discrete  values  as  a  function  of  depth  (with  linear 
interpolation  between  the  given  depths)  or  deter¬ 
mined  from  formulas  as  in  problem  3,  below.  Sepa¬ 
rate  scattering  phase  functions  are  used  for  the 
particles  and  for  the  water  itself.  Variants  of  this 
code  have  been  used  for  a  number  of  studies  of 
radiative  transfer  in  the  ocean.12-17 

The  sea-surface  roughness  is  modeled  using  the 
Cox  and  Munk18  surface  slope  distribution  for  a  given 
wind  speed.  The  effect  of  the  surface  roughness  is 
not  simulated  exactly  because  the  possibility  of  shad¬ 
owing  of  one  facet  by  another  is  ignored.  Multiple 
scattering,  however,  is  included:  e.g. ,  if  a  downward- 
moving  photon  in  the  atmosphere  encounters  the  sea 
surface  and  is  still  moving  downward  after  reflection, 
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it  will  undergo  a  second  interaction  with  the  sea 
surface.  One  important  aspect  of  this  model  is  the 
proper  use  of  photon  weights  to  account  for  the  fact 
that  not  all  facets  are  oriented  in  such  a  manner  .is  to 
be  able  to  interact  with  an  incident  photon,  i.e.,  facets 
with  normals  making  an  angle  less  than  90°  to  the 
direction  of  the  incident  photon.  The  sequence  of 
events  during  an  interaction  with  the  surface  follows. 
From  Cox  and  Munk,  the  probability  that  the  x  and  y 
components  of  the  surface  slope,  zx  and  zy,  respec¬ 
tively,  are  within  zx  ±  '/2d zx  and  zY  ±  Vid zy  is 


p(zx,  zy)dzxdzy 


or 


p{0n,  <J>„  )d0„d<t>„ 


1  /  tan2  <}>„ 

— 2  exp  - - — 

mrt  \  o'2 


x  tan  <J>„  sec2  4>„d0nd<J)n . 


where 


<r2  =  0.003  +  0.0051217. 


Here  U  is  the  wind  speed  in  meters  per  second,  4>n  is 
the  angle  between  the  normal  to  the  facet  and  the 
normal  to  the  level  surface,  and  0„  is  the  azimuth  of 
the  normal.  Given  random  numbers  p8j)  and  p*  on 
the  unit  interval  (0,  1),  the  model  finds  0„  and  4>„  from 


0„  =  2irpen, 


1  / 

wj0  eM- 


tan2  rf>„' 


|tan  <)>„’  sec2  <J>„'d<J),'. 


The  photon  interacting  with  the  surface  is  given  the 
weight 


W  = 


cos  0)  sec  4>„ 


p(zx,  zy)cos  w  sec  <bndzxdzy 

cos<d>0 


and  r  +  dr;  v  =  3  is  used  in  the  computations.  The 
aerosol  total-scattering  coefficient  at  each  altitude  is 
proportional  to  k~p,  where  P  =  u  -  2;  however,  P  = 
0.75  fits  Elterman’s  data  better.  When  a  photon 
interacts  with  the  atmosphere,  the  scattering  angle  is 
chosen  from  either  the  molecular  or  aerosol  phase 
functions  based  on  the  ratio  of  their  scattering  coeffi¬ 
cients  for  the  layer  in  which  the  interaction  takes 
place. 

When  inelastic  processes  are  to  be  included,  the 
above  code  is  operated  at  the  excitation  wavelength,^, 
to  determine  the  excitation  radiance  distribution. 
This  is  used  as  input  to  a  second  Monte  Carlo  code 
that  computes  the  light  field  at  the  wavelength  of 
interest.17  As  with  the  elastically  scattered  radia¬ 
tion,  the  goal  is  to  determine  the  irradiances  of  the 
inelastically  scattered  radiation.  This  is  a  consider¬ 
able  simplification  because  the  solution  can  be  ef¬ 
fected  by  working  with  the  azimuthally  averaged 
radiance  at  A;  i.e.,  only  the  azimuthally  averaged 
radiative  transfer  equation  need  be  solved.  The 
details  of  this  formulation  are  given  in  Appendix  A. 

D.  Model  MC2  (Monte  Carlo  2,  Kattawar) 

This  model  also  simulates  a  coupled  ocean-atmo¬ 
sphere  system.  The  Monte  Carlo  code  relies  heavily 
on  several  variance-reducing  schemes  to  increase 
computational  efficiency.  We  give  only  a  brief  de¬ 
scription  of  one  of  the  most  useful  ones.  The  use  of 
statistical  weights  allows  us  to  treat  each  photon 
history  as  a  packet  of  photons  rather  than  as  a  single 
photon.  Photons  are  never  allowed  to  escape  from 
the  ocean-atmosphere  system.  The  method  of  forced 
collisions  is  used,  whereby  we  sample  from  a  biased 
distribution  that  ensures  a  collision  along  the  path, 
and  the  weight  is  then  adjusted  appropriately  to 
unbias  the  result.  The  way  this  is  done  is  as  follows. 
Suppose  one  wants  to  compute  the  expectation  value 
{ f )  of  some  function  f  oi  a  random  variable  x,  using  a 
probability  density  function p(x).  By  definition, 

(f)  =  J  f{x)p{x) dx. 


where  <d  is  the  angle  of  incidence  upon  the  chosen 
facet.  The  weight,  W,  accounts  for  sampling  from 
p(zx,  zy)  even  though  all  facets  are  not  visible  to  the 
photon. 

The  atmospheric  part  of  the  model  consists  of  fifty 
1-km  layers  with  both  molecular  and  aerosol  scatter¬ 
ing.  The  vertical  distribution  of  the  optical  proper¬ 
ties  is  taken  from  Elterman.19  The  aerosol  phase 
function  at  the  given  wavelength  is  determined  from 
Mie  theory20  with  Deirmendjian’s  Haze  C  size  distri¬ 
bution21 

dn(r)  1 
— ; —  *  — - 1 
dr  r1"1-1 

where  r  is  the  particle  radius  and  d n(r)  is  the  number 
of  particles  per  unit  volume  with  radius  between  r 


However,  if  we  want  to  sample  from  the  density 
function,  p(x),  then 

if)  =  J  f(x)j^7jp(x)  dx  =  J  f(x)w(x)p(x)  dx, 

where  w{x)  =  p(x)/p(x)  is  called  the  statistical  weight. 
The  variance  a2  of  f(x)w(x)  when  sampling  from  the 
biased  distribution  is  given  by 

<r2[  f{x)w(x)}  =  J  [  f(x)w(x)  -  (f)fp{x) dx. 

Although  this  method  appears  straightforward,  it 
does  have  pitfalls.  If  the  weight  can  have  values  that 
exceed  unity,  then  one  can  have  a  variance  that  far 
exceeds  the  variance  in  the  unbiased  sampling. 
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Therefore  extreme  caution  must  be  used  when  apply¬ 
ing  this  method.  It  should  be  noted  that  this  is  a 
very  powerful  method  for  studying  perturbation  ef¬ 
fects,  because  several  processes  can  be  simulta¬ 
neously  emulated  with  the  same  set  of  photon  histo¬ 
ries. 

Now  consider  the  technique  of  forced  collisions,  in 
which  photons  are  never  allowed  to  escape  the  me¬ 
dium.  Let  t*  denote  the  optical  path  length  to  a 
boundary.  To  ensure  that  the  photon  never  escapes, 
we  sample  the  path  length  according  to  the  probabil¬ 
ity  density  function 


p(T)dt  = 


exp(-T)dT 
1  -  exp(-Th)  ’ 


0  <  T  <  Tj,. 


The  weight  now  has  to  be  multiplied  by  [  1  -  exp(  -  t4  )  j 
to  remove  the  bias.  It  should  be  noticed  that  this 
factor  is  always  less  than  unity  and  should  produce  a 
smaller  variance  than  that  produced  when  using 
unforced  sampling.  Histories  are  terminated  only 
when  the  statistical  weight  falls  below  some  specified 
value. 

When  an  interaction  occurs,  the  packet  weight  is 
multiplied  by  the  single-scattering  albedo,  u>0,  which 
gives  the  fraction  of  photons  that  can  continue  to 
scatter.  The  level  air-water  interface  is  modeled  by 
using  the  appropriate  Fresnel  reflection  and  transmis¬ 
sion  coefficients.  A  random  number  is  chosen  at  this 
stage  to  determine  whether  the  photon  is  transmitted 
or  reflected. 

Radiances  are  obtained  over  detectors  that  have 
finite  solid  angles.  However,  statistical  estimation 
can  be  used  to  give  true  continuum  radiance  values 
where  no  directional  averaging  is  done.  This  model 
can  simulate  inelastic  scattering;  the  details  are  given 
in  Kattawar  and  Xu.22  The  Monte  Carlo  method  has 
also  been  extended  to  include  the  full  Stokes  vector 
treatment  of  polarization;23-26  these  papers  show  that 
substantial  errors  can  occur  if  polarization  is  ne¬ 
glected. 

E.  Model  MC3  [Monte  Carlo  3  (Morel  and  Gentili)] 

This  Monte  Carlo  model  is  similar  to  those  described 
by  Plass  and  Kattawar2728  and  by  Gordon  and 
Brown.29  It  is  designed  to  simulate  the  radiance 
distribution  at  any  level  in  the  atmosphere  and  in  the 
ocean.  Between  these  two  media,  a  wind-roughened 
interface  is  modeled  with  the  isotropic  Gaussian 
distribution  of  sea-surface  slopes,  as  discussed  under 
model  MCI.  The  probability  of  occurrence  of  the 
various  slopes  is  modified  when  considering  nonverti- 
cally  incident  photons.  This  photon-facet  interac¬ 
tion  is  modeled  as  in  Plass  et  al.;30  it  does  not  account 
for  the  possible  occultation  of  a  facet  by  an  adjacent 
one.  Transmitted  and  reflected  photon  packets  re¬ 
sulting  from  interaction  with  the  air-water  surface 
are  weighted  according  to  Fresnel’s  law  (including  the 
possibility  of  total  internal  reflection).  According  to 
the  problem  under  investigation,  photon  packets  are 
introduced  at  the  top  of  the  atmosphere,  or  just  above 
(or  below)  the  ocean  surface.  For  specific  problems 


involving  deep  levels,  packets  can  be  reintroduced  at 
intermediate  depths  inside  the  water  body,  according 
to  a  directional  distribution  that  reproduces  the 
downward  radiance  field  as  resulting  from  a  previous 
Monte  Carlo  run.  The  bottom  boundary  is  either  an 
infinitely  thick  absorbing  layer,  in  which  photons  are 
lost  from  the  system,  or  a  Lambertian  reflecting 
bottom  of  a  given  albedo,  from  which  weighted  pho¬ 
ton  packets  are  reflected. 

After  each  collision,  the  weight  of  each  photon 
packet  is  multiplied  by  the  local  value  of  oj0  that  is 
pertinent  to  the  altitude  or  the  depth,  to  account  for 
its  partial  absorption.  A  packet  history  is  termi¬ 
nated  when  its  weight  falls  below  a  predetermined 
value,  typically  1  x  10  ~6.  For  each  collision  a  ran¬ 
dom  number  on  the  unit  interval  is  compared  with 
the  local  value  of  the  ratio  of  the  molecular  scattering 
coefficient  to  the  total  scattering  coefficient  to  deter¬ 
mine  if  the  scattering  event  will  be  of  molecular  type 
(air  or  water  molecules)  or  is  due  to  an  aerosol  or 
hydrosol  particle.  The  appropriate  phase  function  is 
then  used  to  determine  the  scattering  angle;  the 
orientation  of  the  scattering  plane  is  chosen  at  ran¬ 
dom  on  the  interval  (0,  2-rr).  The  number  of  photons 
initiated  depends  on  the  single-scattering  albedo  value, 
so  as  to  control  the  stochastic  noise  in  the  computed 
radiometric  quantities  (details  can  be  found  in  Morel 
and  Gentili31-32).  The  model  is  operated  for  its  oce¬ 
anic  segment  with  the  optical  properties  as  specified 
in  Section  3.  For  the  atmospheric  segment,  fifty 
1-km-thick  layers  are  considered,  with  specified  val¬ 
ues  for  Rayleigh  and  aerosol  scattering  and  for  ozone 
absorption  as  in  Elterman.19  The  aerosol  phase 
function  (as  computed  by  Mie  scattering  theory)  for 
the  maritime  aerosol  model  defined  by  the  Radiation 
Commission  of  the  International  Association  of  Meteo¬ 
rology  and  Atmospheric  Physics  is  used;  see  the 
models  of  Tanre  et  al.33  and  Baker  and  Frouin.34 


F.  Model  MC4  [Monte  Carlo  4  (Reinersman)) 

This  model  is  intended  primarily  for  simulation  of  the 
radiance  distribution  above  and  just  below  the  sur¬ 
face,  and  for  simulation  of  irradiances  with  the  first 
five  mean  free  paths  of  the  surface.  The  model  is 
based  on  techniques  described  by  Kirk.35  The  model 
atmosphere  is  composed  of  50  layers,  each  character¬ 
ized  by  separate  Rayleigh  and  particulate  scattering 
coefficients  and  an  albedo  of  single  scattering,  as 
given  by  Elterman.19  Weighted  photon  beams  are 
projected  into  the  atmosphere  from  the  atmosphere- 
space  boundary,  and  a  colllision  is  forced  somewhere 
in  the  atmosphere  along  this  original  trajectory. 
The  attenuated  beam,  which  is  the  weight  of  the 
original  beam  less  the  portion  lost  to  scattering  and 
absorption,  strikes  the  sea  surface  at  the  angle  of  the 
original  trajectory.  Beam  losses  that  are  due  to 
absorption  and  scattering  take  place  at  the  point  of 
collision.  There  the  absorbed  portion  is  lost  and  the 
scattered  portion  exits  the  collision  point  in  another 
single,  weighted  beam.  A  random  number  is  com¬ 
pared  with  the  ratio  of  the  Rayleigh  scattering  cross 
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section  to  the  total  scattering  cross  section  to  deter¬ 
mine  the  type  of  volume  scattering  function  govern¬ 
ing  the  scattering  event.  In  the  case  of  an  aerosol 
scattering,  a  two-term  Henyey-Greenstein  phase  func¬ 
tion  is  used  to  determine  the  scattering  angle.36 
Otherwise  the  angle  is  determined  by  a  Rayleigh 
phase  function.37  Once  the  trajectory  of  the  scat¬ 
tered  portion  of  the  beam  is  calculated,  the  distance 
from  the  point  of  collision  to  the  next  encountered 
interface  (air-water  or  air-space)  is  determined. 
A  new  collision  is  forced  somewhere  along  this  trajec¬ 
tory,  and  the  process  is  repeated  until  the  '  eight  of 
the  scattered  portion  of  the  beam  falls  below  a  preset 
minimum  fraction  of  the  original  beam  weight.  This 
minimum  traceable  weight  is  set  to  1  x  10  6  of  the 
original  beam  weight  for  the  simulations  presented 
below. 

Some  of  the  scattered  trajectories  encounter  the 
atmosphere-space  boundary  and  are  forgotten;  the 
others  impinge  on  the  sea  surface.  For  the  latter, 
the  angle  of  incidence  depends  on  the  nadir  angle  of 
the  ray  and  the  slope  of  the  sea  surface.  The  direc¬ 
tions  of  the  reflected  and  refracted  rays  are  deter¬ 
mined  geometrically,  and  the  weights  of  the  rays  are 
calculated  from  the  Fresnel  formula.  Although  wave 
shadowing  is  neglected,  multiple  surface  interactions 
may  occur.  A  reflected  ray  that  is  still  projected 
downward,  or  a  transmitted  ray  that  is  still  projected 
upward,  must  encounter  the  sea  surface  again  imme¬ 
diately,  without  an  intervening  trajectory.  Ray  tra¬ 
jectories  resulting  from  reflection  are  followed  in  the 
original  manner.  Transmitted  portions  of  the  beams 
are  followed  similarly  until  encountering  the  bottom 
or  the  sea  surface,  or  until  they  are  diminished  to  less 
than  the  minimum  traceable  weight.  Those  beams 
striking  the  bottom  are  lost;  those  beams  that  are 
incident  upon  the  sea  surface  from  below  are  again 
subjected  tc  the  reflection  and  transmission  calcula¬ 
tions. 

G.  Model  MC5  [Monte  Carlo  5  (Stavn)J 
The  Naval  Research  Laboratory  optical  model  (re¬ 
ferred  to  as  the  NORDA  or  NOARL  optical  model  in 
earlier  publications)  uses  standard  Monte  Carlo  tech¬ 
niques.1328’35  At  each  scattering  event,  a  random 
number  is  used  to  determine  if  the  scattering  is  due  to 
molecular  water,  quartzlike  particulates,  algae,  or 
organic  detritus;  the  volume  scattering  functions  of 
these  components  are  treated  separately,  rather  than 
using  an  average  volume  scattering  function.  The 
model  includes  the  effects  of  Raman  scattering.  If  a 
photon  collision  results  in  inelastic  scattering  (as 
determined  by  comparing  a  random  number  to  the 
appropriate  optical  properties  of  the  medium),  the 
wavelength  is  shifted  by  an  amount  corresponding  to 
the  mean  wave-number  shift  of  3357  cm-1,  corre¬ 
sponding  to  Raman  scatter  by  water  molecules.  The 
finite  bandwidth  of  the  Raman-shifted  light  is  taken 
into  account  by  averaging  over  10-nm  bandwidths 
(roughly  corresponding  to  current  oceanographic  in¬ 
struments);  details  of  this  averaging  are  described  in 


Stavn  and  Weidemann.38  39  For  the  simulation  of 
problem  7,  below,  it  was  assumed  that  the  Raman 
scattering  occurs  in  a  very  narrow  waveband.  The 
photons  are  tallied  into  zonal  bands,  as  is  convenient 
for  computation  of  irradiances  and  the  nadir-viewing 
radiance. 

There  is  no  atmosphere  per  se  implemented  in  the 
model.  Atmospheric  transmittances  of  solar  irradi- 
ance  needed  for  simulations  are  obtained  from  the 
nonlayered  atmospheric  model  of  Brine  and  Iqbal.40 
The  model  determines  the  skylight  radiance  pattern 
from  the  empirical  model  of  Harrison  and  Coombes.5 
The  present  version  of  the  code  handles  only  homoge¬ 
neous  waters. 

3.  Canonical  Problems 

We  now  define  several  canonical,  or  standard,  prob¬ 
lems  for  solution  by  underwater  radiative  transfer 
models.  Models  claiming  to  provide  realistic  simula¬ 
tions  of  the  oceanic  optical  environment  should  be 
able  to  solve  these  problems  and  provide  output  that 
is  at  least  as  accurate  as  the  data  obtainable  by 
presently  available  instrumentation.  In  brief,  these 
problems  are 

Problem  1:  An  unrealistically  simple  problem. 

Problem  2:  A  base  problem  using  realistic  inher¬ 
ent  optical  properties  for  the  sea  wa¬ 
ter. 

Problem  3;  The  base  problem  but  with  stratified 
water. 

Problem  4:  The  base  problem  but  with  atmo¬ 
spheric  effects. 

Problem  5:  The  base  problem  but  with  a  wind¬ 
blown  sea  surface. 

Problem  6:  The  base  problem  but  with  a  finite- 
depth  bottom. 

Problem  7:  A  problem  involving  Raman  scatter¬ 
ing. 

In  each  of  these  problems,  the  water  body  is  taken 
to  be  horizontally  homogeneous.  The  real  index  of 
refraction  of  the  water  is  n  =  1.340.  The  depth 
below  the  surface  can  be  specified  by  either  the 
nondimensional  optical  depth  t  or  by  the  geometric 
depth  z  in  meters.  The  base  problem  2  assumes  that 
(a)  the  air-water  surface  is  flat;  (b)  the  water  is 
homogeneous  and  infinitely  deep;  (c)  there  is  no 
atmosphere,  i.e.,  the  sky  is  black;  (d)  the  sun  is  a  point 
light  source  located  at  a  zenith  angle  of  0sun  =  60°,  (e) 
the  sun  provides  a  spectral  irradiance  just  above  the 
sea  surface  of  magnitude  =  1  W  m  2  nm  1  on  a 
surface  perpendicular  to  the  sun’s  rays  (which  gives 
Ed  =  0.5  W  m~2  nm-1  for  0sun  =  60°);  (f)  there  is  no 
inelastic  scattering  or  other  source  of  light  within  the 
water  body;  (g)  the  angular  scattering  properties  of 
the  water  are  characteristic  of  natural  hydrosols;  and 
(h)  the  water  is  either  highly  scattering  (w0  =  0.9)  or 
highly  absorbing  (w0  =  0.2).  The  other  problems  are 
defined  by  exceptions  to  these  assumptions.  The 
specific  problem  definitions  are  as  follows. 
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Problem  1 .  A  Rayleigh  phase  function 

<i>'  -*  4>)  =  Lw  =  (1  +  cos2  <|r)  (3) 

is  used  to  describe  the  angular  scattering  properties 
of  the  water.  The  scattering  angle,  il»,  is  related  to 
the  incident  (jjl',  <f>')  and  scattered  (jx,  4>)  directions  by 

rj>  =  cos~l[p.p,’  +  (1  -  U2)1/2(l  -  M-'2)1  2COs(<i>  - 

This  phase  function,  which  is  plotted  in  Fig.  1(b),  is 
similar  to  that  of  pure  sea  water.  The  Rayleigh 
phase  function  is  a  well-behaved  function  of  the 
scattering  angle,  <Ji,  and  presents  no  numerical  difficul¬ 
ties  in  its  treatment;  we  therefore  consider  this  an 
easy  problem  for  numerical  modeling.  Note  that  (3^ 
satisfies  the  normalization 

2ir  I  p(tjr)sin  rjrdrjr  =1.  (4) 

Jo 

Both  highly  scattering  (co0  =  0.9)  and  highly  absorb- 


(b) 

Fig.  1.  (a)  Inherent  optical  properties  as  a  function  of  depth  for 
problem  3.  Coefficients  a,  b,  and  c  have  units  of  inverse  meters;  ioo 
is  dimensionless;  (b)  scattering-phase  function  for  pure  sea  water, 
bu,;  for  particles,  fjp;  and  for  problem  3  at  depths  of  2  =  0,  17,  and 
60  m. 


ing  (u>o  =  0.2)  cases  are  considered  for  the  Rayleigh 
phase  function. 

Problem  2.  This  base  problem  uses  a  phase  func¬ 
tion  that  is  typical  of  oceanic  waters.  The  total 
volume  scattering  function  (VSF)  p  is 

P  =  Pa  +  Pp, 

where  subscripts  w  and  p  refer  to  pure  sea  water  and 
to  particles,  respectively.  The  total  phase  function  (3 
therefore  can  be  expressed  as 

bw  -  bD  - 

P  =  fPa  +  fPp-  (5) 

This  total  (3  must  satisfy  the  normalization  (4),  which 
is  the  case  if  (3a,  and  |3P  are  each  normalized. 

The  particle  phase  function,  Pp,  is  defined  from 
three  VSF’s  measured  by  Petzold41  in  San  Diego 
harbor.  The  VSF  for  pure  sea  water42  was  first 
subtracted  to  find  the  three  particle  VSF’s.  Then 
the  scattering  coefficient  of  pure  sea  water43  (bw  = 
0.00231  m_1  at  \  =  530  nm,  the  wavelength  of 
Petzold ’s  data)  was  subtracted  from  the  respective 
scattering  coefficients  computed  by  Petzold  (b  = 
1.205,  1.536,  and  1.824  m_1  for  the  three  VSF’s)  to 
find  the  particle-scattering  coefficient,  bp,  for  each 
VSF.  The  three  particle  phase  functions  were  then 
computed  with  these  bp's,  and  the  mean  value  of  the 
three  Pp’s  was  computed  at  each  scattering  angle. 
This  mean  pp[ijr)  becomes  infinite  at  i{j  =  0,  if  it  is 
assumed  that  Pp(»|j)  ~  i) rm  as  i|»  -*  0,  where  m  =  1.346 
is  the  negative  of  the  slope  of  log  |3P(»|/)  versus  log  r|»  at 
the  two  smallest  tabulated  scattering  angles  (4»  = 
0.10°  and  0. 12589°).  When  this  functional  form  of  (3P 
was  used  to  integrate  analytically  2irpp(r)r)sin  r]r  from 
t|/  =  0  to  =  0.10°,  and  the  trapezoidal  rule  was  used 
to  integrate  from  «|»  =  0. 10°  to  i|/  =  180°,  the  normaliza¬ 
tion  integral  (4)  gave  the  value  1.006449.  We  thus 
divided  the  mean  pp  by  1.006449  to  obtain  the  values 
shown  in  Table  2.  The  particle  phase  function  Pp(i|i) 
is  then  defined  to  be  the  tabulated  values,  with  linear 
interpolation  to  be  used  between  the  tabulated  values 
and  with  pp(»fi)  ■  pp(0.12589°J  (0.12589°/ v)))1 346  for 
ip  <  0.12589°.  The  resulting  Pp(i|/)  is  defined  for  all 
and  exactly  satisfies  the  normalization  condition  (4). 
This  Pp  is  plotted  in  Fig.  1(b). 

Moreover,  because  bw  =  0.00231  mr1  is  much  less 
than  bp(  >  1.2  m'1  for  each  of  the  Petzold  VSF’s),  it  is 
reasonable  to  neglect  the  contribution  of  the  water, 
Pu,,  to  the  total  phase  function  of  Eq.  (5).  This 
omission  creates  an  error  of  at  most  a  few  percent  in  p 
even  at  backscattered  directions  (i(>  >  90°).  We  there¬ 
fore  define  the  total  phase  function  for  problem  2  to 
be  just  the  particle  phase  function  as  defined  above; 
P(d>)  =  (3p(»li).  This  p  is  representative  of  phase 
functions  measured  in  ocean  waters  with  typical 
particle  concentrations  and,  because  of  its  highly 
peaked  behavior  at  small  i|«,  can  be  expected  to  test  the 
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Tabta  2.  Ptiaae  Function  Values  Used  in  (Mining  the  Particulate  Phase 
Function 


Scattering 

Angle 

(deg) 

Phase 

Function" 

(sr-‘l 

Scattering 

Angle 

(deg) 

Phase 

Function 

(sr_1) 

0.10000 

1.76661+3 

50.0 

2.27533-2 

0.12589 

1.29564+3 

55.0 

1.69904-2 

0.15849 

9.50172+2 

60.0 

1.31254  -2 

0.19953 

6.99092+2 

65.0 

1.04625-2 

0.25119 

5.13687+2 

70.0 

8.48826-3 

0.31623 

3.76373+2 

75.0 

6.97601-3 

0.39811 

2.76318+2 

80.0 

5.84232-3 

0.50119 

2.18839+2 

85.0 

4.95306-3 

0.63096 

1.44369+2 

90.0 

4.29232-3 

0.79433 

1.02241+2 

95.0 

3.78161-3 

1.0000 

7.16082  +  1 

100.0 

3.40405-3 

1.2589 

4.95803+1 

105.0 

3.11591-3 

1.5849 

3.39511  +  1 

110.0 

2.91222-3 

1.9953 

2.28129+1 

115.0 

2.79696  -3 

2.5119 

1.51622  +  1 

120.0 

2.68568  -3 

3.1623 

1.00154+1 

125.0 

2.57142-3 

3.9811 

6.57957 

130.0 

2.47603-3 

5.0119 

4.29530 

135.0 

2.37667-3 

6.3096 

2.80690 

140.0 

2.32898-3 

7.9433 

1.81927 

145.0 

2.31308-3 

10.0 

1.15257 

150.0 

2.36475-3 

15.0 

4.89344-1 

155.0 

2.50584-3 

20.0 

2.44424-1 

160.0 

2.66183-3 

25.0 

1.47151-1 

165.0 

2.83472-3 

30.0 

8.60848-2 

170.0 

3.03046-3 

35.0 

5.93075-2 

175.0 

3.09206-3 

40.0 

4.20985-2 

180.0 

3.15366-3 

45.0 

3.06722-2 

“The  notation  n  ±  e  =  n  x  10-*. 


numerical  models’  abilities  to  handle  realistic  phase 
functions.  Both  highly  scattering  and  highly  absorb¬ 
ing  cases  are  considered  for  this  phase  function. 

Problem  3.  This  problem  is  designed  to  test  the 
models’  abilities  to  compute  light  fields  in  highly 
stratified  water.  The  water  stratification  is  specified 
as  follows.  The  particulate  absorption  and  scatter¬ 
ing  coefficients  are  taken  to  be 

ap  =  0.04C0602,  (6a) 

bp  =  0.33C0620,  (6b) 

respectively,  where  C  is  the  chlorophyll  (pigment) 

concentration.  When  C  is  in  milligrams  per  inverse 

meters  cubed,  ap  and  bp  are  in  inverse  meters.  The 
absorption  representation  (6a)  is  based  on  Prieur  and 
Sathyendranath44  at  a  wavelength  of  X.  =  500  nm. 
The  scattering  representation  (6b)  is  based  on  Gordon 
and  Morel45  with  X  =  500  nm  and  assuming  that 
bp(  X )  ~  X  - 1 .  The  pigment  profile  with  depth  is  based 
on  Lewis  et  al.46  and  consists  of  a  Gaussian  plus  a 
constant  background: 


C{z)  =  C0  + 


(7a) 


Platt  and  Sathyendranath47  show  that  Eq.  (7a)  with 


the  parameter  values 


C0  =  0.2  mg  m"3,  (7b) 

s  =  9  m,  (7c) 

2mai  =  17  m,  (7d) 

h  =  144  mg  m'2,  (7e) 


fits  data  from  the  Celtic  Sea  in  May  very  well.  We 
therefore  adopt  Eq.  (7)  as  a  reasonable  model  for  C(r). 
When  Eq.  (7)  is  used  in  Eq.  (6),  the  particulate 
absorption  and  scattering  coefficients,  and  hence  all 
inherent  optical  properties,  become  functions  of  depth. 
The  absorption  and  scattering  coefficients  for  pure 
sea  water  at  X  =  500  nm  are  given  by43 

aw  =  0.0257  m'1  (8a) 

bw  -  0.0029  m-‘.  (8b) 

When  the  chlorophyll  concentration  is  low,  scatter¬ 
ing  by  pure  sea  water  makes  a  significant  contribu¬ 
tion  to  the  total  scattering  at  large  scattering  angles 
(almost  1/2  when  C  =  C0  and  ij»  =  180°).  Therefore, 
for  this  problem  it  is  necessary  to  use  Eq.  (5)  to 
determine  the  total  phase  function  from  the  phase 
functions  for  pure  sea  water,  (3^,  and  for  particles,  [}p, 
as  were  defined  in  problems  1  and  2.  The  phase 
function  is  now  a  function  of  depth,  as  is  the  scatter- 
ing-to-attenuation  ratio 

=  -  bw  +  bP^ 

°>0~  c  aw  +  ap(z)  +  bw  +  bp(z ) 

Figure  1(a)  shows  a,  b,  c,  and  u>o  as  functions  of  depth 
for  problem  3,  and  Fig.  1(b)  shows  the  phase  func¬ 
tions  at  selected  depths. 

Problem  4.  This  problem  is  the  same  as  problem  2 
with  u)q  =  0.9,  except  that  atmospheric  effects  are 
included.  The  sky  is  no  longer  black  but  rather  has  a 
radiance  distribution  that  describes  the  atmosphere’s 
scattering  and  absorption  effects  on  sunlight.  The 
incident  solar  irradiance,  E±  =  1 W  m~2  nm-1,  is  now 
applied  at  the  top  of  the  atmosphere.  The  atmo¬ 
spheric  optical  effects  are  defined  by  Elterman’s19 
aerosol  and  Rayleigh-scattering  optical  thicknesses  at 
X  =  500  nm: 


Vrosol  =  0.264, 
^Rayleigh  “  0.145. 


Because  the  numerical  models  incorporate  atmo¬ 
spheric  effects  in  various  ways,  a  more  detailed 
specification  of  the  atmosphere  is  not  made. 

Problem  5.  This  problem  is  the  same  as  problem  2 
with  u)0  =  0.9,  except  that  the  effects  of  a  windblown 
sea  surface  are  included.  The  surface  waves  are 
statistically  specified  as  having  a  wave  slope  standard 
deviation  of  cr  =  0.2  in  the  Cox-Munk18  capillary- 
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wave  spectrum 

a2  =  0.003  +  0.00512 U, 

where  U  is  the  wind  speed  in  meters  per  second. 
Thus  cr  =  0.2  corresponds  to  a  wind  speed  of  "  =  7.23 
ms1.  The  solar  zenith  angle  is  taken  to  oe  0BUn  = 
80°. 

Problem  6.  This  problem  is  the  ^ame  as  problem 
2,  except  that  a  finite-depth  bottom  is  imposed.  The 
bottom  is  taken  to  be  an  opaque,  Lambertian  reflect¬ 
ing  surface  at  depth  t  =  5.  This  surface  has  an 
irradiance  reflectarne(Eu  Ed )  of  0.5.  Such  a  surface 
is  a  reasonable  n.odel  of  a  light-colored,  sandy  bot¬ 
tom. 

Problem  7.  This  problem  is  for  use  in  comparing 
models  that  include  the  effects  of  Raman  scattering 
by  water  molecules.  The  wavelength  of  excitation  is 
taken  to  be  \ex  =  417  nm,  and  all  light  that  is  Raman 
scattered  at  417  nm  is  assumed  to  shift  to  K  =  486 
nm.  The  Rayleigh  phase  function,  Eq.  (3),  is  used 
for  elastic  scattering.  The  phase  function  for  Ra¬ 
man  scattering  is48 

3  1  +  3p  /  1-p  \ 

PiuJ'l')  -  ^  i  +  2p  I1  +  i  +  3p  cos  *)  ’  !9) 

where  p  is  the  depolarization  ratio.  For  this  prob¬ 
lem,  we  use  p  =  0.17  and  take  the  total  Raman 
scattering  coefficient  bf^m  equal  to  the  elastic-scatter¬ 
ing  coefficient  of  the  water  itself,  i.e.,  bf^m  =  bu..  The 
absorption  and  elastic-scattering  coefficients  of  pure 
sea  water  at  the  wavelengths  in  question  as  taken 
from  Smith  and  Baker43  are 

<U417)  -  0.0156  m  », 
bj 417)  =  0.0063  m1, 
aJ486)  =  0.0188  m1, 

6,„(486)  =  0.0032  m1. 

Considering  the  way  in  which  Smith  and  Baker 


inferred  au  from  irradiance  data,  it  is  assumed  that 
^Ram  is  already  included  in  the  value  of  au .  Thus  the 
total  beam  attenuation  coefficient  at  each  wavelength 
is  just  aw  +  bu.  A  unit  irradiance  Ek  is  incident  at 
the  sea  surface  upon  a  plane  normal  to  the  solar  beam 
at  the  excitation  wavelength  Aei  =  417  nm.  There  is 
no  atmosphere  and  no  solar  irradiance  is  incident 
upon  the  sea  surface  at  X  =  486.  The  resulting 
irradiances  at  486  nm  are  those  that  would  be  solely 
because  of  inelastic  scattering  from  417  nm.  The 
solar  zenith  angle  is  60°  and  the  air-water  surface  is 
flat. 

Table  3  summarizes  the  various  canonical  prob¬ 
lems. 

4.  Model  Comparison* 

Although  the  models  generally  compute  the  radiance 
L,  the  quantities  most  often  used  in  oceanic  optics  are 
various  irradiances.  These  irradiances  are  defined 
by  weighted  integrations  of  the  radiance  distribution 
over  the  upward  and  downward  hemispheres  of  direc¬ 
tions,  as  shown  in  Table  1,  and  are  easily  obtained 
from  computed  radiances.  The  nadir-viewing  radi¬ 
ance,  Lu,  is  the  radiance  seen  by  a  sensor  pointed 
straight  down  (in  the  nadir  direction);  L„  is  important 
in  remote-sensing  studies.  The  ability  of  a  numeri¬ 
cal  model  to  accurately  compute  the  irradiances  and 
nadir  radiance  is  a  measure  of  its  utility  for  many 
oceanographic  studies. 

Models  II  and  DO  compute  all  quantities  with  equal 
accuracy.  However,  the  Monte  Carlo  models  MC1- 
MC5  compute  upwelling  quantities  (e.g.,  Eu,  Eou,  or 
Lu  )  with  less  accuracy  than  downwelling  quantities 
(e.g.,  Erf  or  E^).  This  is  because  most  of  the  simu¬ 
lated  photons,  all  of  which  are  initially  heading 
downward,  continue  to  head  downward  and  thereby 
contribute  to  Ed  or  E<ld .  However,  only  the  relatively 
few  photons  that  are  scattered  into  upward  directions 
can  contribute  to  Eu,  Eou,  or  Lu;  fewer  photons  means 
greater  statistical  fluctuations  in  the  computed  val¬ 
ues. 

Also,  for  a  given  initial  number  of  photons,  the 


Tabt*  3.  Summary  of  tha  Canonical  Problem* 


Problem 

Parameter 

1 

Easy 

Problem 

2 

Base 

Problem 

3 

Stratified 

Water 

4 

Atmospheric 

Effects 

5 

Windblown 

Surface 

6 

Bottom 

Effects 

7 

Raman 

Scattering 

Albedo,  too 

0.9, 0.2 

0  9,  0.2 

Depth 

0.9 

0.9 

02 

0  29  at  417  nm 

Phase 

Rayleigh 

Particle 

dependent 

Depth 

Particle 

Particle 

Particle 

0. 15  at  486  nm 

Eqs.  (3)  and  (9) 

function 

Eq.  (3) 

Table  2 

dependent 

Table  2 

Table  2 

Table  2 

Air-water 

Flat 

Flat 

Flat 

Flat 

Capillary 

Flat 

Flat 

surface 
Diffuse  sky 

0 

0 

0 

Various 

waves 

0 

0 

0 

radiance 

Internal 

0 

0 

0 

models 

0 

0 

0 

Various  models 

sources 

Bottom 

Infinitely 

Infinitely 

Infinitely 

Infinitely 

Infinitely 

Lambertian 

Infinitely 

boundary 

deep 

deep 

deep 

deep 

deep 

at  t  =  5 

deep 
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Monte  Carlo  models  must  settle  for  less  accuracy  at  a 
given  optical  depth  t  in  highly  absorbing  waters 
(small  wo)  than  in  highly  scattering  waters  (large  w()). 
This  is  because  photons  absorbed  before  they  reach 
depth  t  are  not  available  to  be  tallied  in  the  computa¬ 
tion  of  the  radiance  or  irradiance,  whereas  scattered 
photons  can  eventually  reach  depth  t  and  be  tallied. 
In  practice,  the  accuracy  of  the  Monte  Carlo  models  is 
strongly  dependent  on  the  number  of  photon  colli¬ 
sions;  thus  more  photons  must  be  processed  when  wo 
is  small  to  achieve  satisfactory  accuracy.  The  accu¬ 
racy  of  models  II  and  DO  is  independent  of  w0. 

With  the  above  comments  in  mind,  we  selected  Ed , 
Eou,  and  Lu  for  comparison  just  above  the  sea  surface 
and  at  t  =  1,5,  and  10.  Problems  1  and  2  have  both 
higly  scattering  (w0  =  0.9)  and  highly  absorbing 
(wo  =  0.2)  waters. 

Although  it  is  not  possible  to  compare  the  computa¬ 


tional  efficiencies  of  the  various  models  because  they 
were  run  on  a  variety  of  computers,  with  differing 
numbers  of  photons  traced  in  the  Monte  Carlo  codes, 
Table  4  shows  some  representative  execution  times. 
It  should  be  noted  that  the  long  execution  times 
shown  for  some  of  the  Monte  Carlo  codes  are  the 
times  required  for  accurate  radiance  simulations  at 
large  depths.  If  only  irradiances  or  near-surface 
radiances  are  required  for  a  particular  study,  these 
models  can  be  run  for  much  shorter  times.  For 
example,  in  the  simulation  of  problem  3,  output  from 
model  MCI  was  compared  for  run  times  of  180  s  and 
7200  s.  The  Ed  values  throughout  the  euphotic  zone 
(roughly  the  upper  21  m),  as  acc  .mulated  after  180  s, 
were  within  1.5%  of  the  values  obtained  after  7200  s. 
After  180  s,  the  Enu  and  Lu  values  just  below  the 
surface  (at  z  =  0)  were  within  1%  of  their  final  values. 
Deeper  within  the  euphotic  zone,  Eou  and  L„  differed 


Tabto  4.  Raprasantattv*  Execution  Timex,  and  Number*  of  Simulated  Photon*  for  Model*  MC1-MC5 


Problem 

Execution 

Time 

(8) 

Number  of 

Photons 

Initiated 

Number  of 

Photon 

Collisons 

Model  II  (Computer:  Sun  SPARCstation  2,  no  code  optimization) 

1,  a»o  —  0.9 

349  for  t  =  10;  730  for  t  =  20 

l.iuo  =  0.2 

350  for  t  =  10;  733  for  r  =  20 

2,  u>o  ~  0.9 

306  for  t  =  10;  496  for  t  =  20 

2,  uh)  —  0.2 

386  for  t=  10;  711  for  t  =  20 

3 

1 180  for  z  =  60  m 

Model  DO  (Computer. 

Decstation  5 000  240,  no  code  optimization  ) 

1 ,  too  ~  0  9 

5  for  irradiances  only,  2  layers 

1,  too  ~  0.2 

5  for  irradiances  only,  2  layers 

2,  u)q  =  0.9 

9  for  irradiances  only.  2  layers;  435  for  radiances.  2  layers 

2,  <*>o  —  0.2 

9  for  irradiances  only,  2  layers 

3 

171  for  irradiances  only,  25  layers 

Mode!  MCI  (Computer: 

Decstation  5000) 

1,  u>o  =  0.9 

7200 

1.25  x  106 

4.98  x  107 

1,  u>o  —  0.2 

7200 

6.63  x  106 

3.99  x  107 

2,  u>o  —  0.9 

7200 

9.66  x  106 

7.18  x  107 

2,  u>o  =  0.2 

7200 

7.17  x  106 

3.77  x  107 

3 

7200 

7.49  x  106 

8.74  x  107 

Model  MC2  (Computer: 

Vax  9000) 

1,  u>o  =  0.9 

5830 

1.0  x  106 

9.47  x  107 

1,  w#  “  0.2 

530 

1.0  x  106 

7.54  x  107 

2,  wo  =  0.9 

4630 

1.0  x  106 

9.72  x  107 

2,  u>o  =  0.2 

410 

1.0  x  106 

7.85  x  107 

Model  MC3  (Computer: 

Hewlett  Packard  9000/730) 

1,  a»o  ~  0.9 

60000 

10.9  x  106 

6.72  x  108 

1,  wo  —  0.2 

74000 

55.7  x  106 

7.07  x  108 

2,  u>o  —  0.9 

45000 

8.7  x  106 

7.30  x  108 

2,  too  =  0.2 

84000 

63.7  x  10« 

12.10  x  108 

3 

56000 

8.9  x  106 

9.02  x  108 

Model  MC4  (Computer: 

Microvax  III) 

1,  too  =  0.9 

15100 

5.0  x  104 

1.66  x  107 

1,  too  —  0.2 

17700 

1.0  x  106 

1.44  x  107 

2,  wo  =  0.9 

9680 

8  0  x  104 

1.24  x  107 

2,  too  ~  0.2 

10000 

1.2  x  106 

1.02  x  107 

3 

24200 

1.0  x  105 

3.06  x  107 

Model  MC5  (Computer: 

Cray  Y-MP,  no  vectorization) 

1,  wo  =  0.9 

1981  for  t  =  20 

1.0  x  107 

1,  wo  =  0.2 

416  for  t  =  10 

1.0  x  107 

2,  u>o  =  0.9 

2300  for  t  =  20 

1.0  x  107 

2,  too  —  0.2 

389  for  t  =  10 

1.0  x  107 
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by  as  much  as  8%  and  20%,  respectively,  for  the  two 
run  times.  At  a  depth  of  z  =  60  m,  the  differences  in 
the  computed  quantities  for  the  two  times  were  3% 
for  Ed ,  19%  for  Eout  and  a  factor  of  six  for  Lu.  Model 
DO  is  much  more  efficient  for  irradiance  and  nadir  (or 
zenith)  radiance  computations,  than  for  full  radiance 
computations,  because  only  the  azimuthally  averaged 
equation  (i.e.,  the  m  =  0  component  of  the  radiance)  is 
required  to  compute  irradiances  and  nadir  or  zenith 
radiances.  Full  off-nadir  or  off-zenith  radiance  com¬ 
putations  require  the  evaluation  of  additional  azi¬ 
muthal  components.  Strongly  anisotropic  scatter¬ 
ing  also  requires  a  large  number  of  streams. 

We  now  briefly  discuss  the  results  of  the  models’ 
simulations  of  problems  1-7. 

Problem  1.  Figure  2(a)  shows  the  computed  Ed, 
Eou,  and  Lu  for  the  Rayleigh  phase  function  of  prob- 


(a) 


spectral  irrodiance  E  (W  m"’  nm"') 
spectral  rodionce  Lu  (W  m"J  sr"'  nm'') 

(b) 


Fig.  2.  (a)  Ed,  Eou,  and  Lu  as  computed  by  the  various  models  for 
problem  1,  coo  =  0.9;  (b)  the  same  quantities  as  computed  for  the 
case  of  coo  “  0.2.  The  dotted  line  represents  the  air-water 
surface.  Results  from  models  II  and  DO  are  plotted  with  solid 
lines;  models  MC1-MC5  are  plotted  with  dashed  lines.  Depth  t  = 
0  is  in  the  water,  just  below  the  surface,  and  in  air  represents  a 
point  just  above  the  surface. 


lem  1  and  o>0  =  0.9.  In  this  and  subsequent  figures, 
we  plot  the  results  from  the  two  analytic  models,  II 
and  DO,  with  solid  lines;  the  Monte  Carlo  results  are 
plotted  with  dashed  lines.  This  makes  it  easy  to  see 
that,  in  most  instances,  the  Monte  Carlo  results  are 
distributed  to  either  side  of  the  analytic  results, 
which  are  usually  indistinguishable  in  the  figures. 

We  first  note  in  Fig.  2(a)  that  all  models  predict 
nearly  the  same  values  for  a  given  quantity,  although 
there  is  a  detectable  spread  in  Lu  values  that  is  due  to 
Monte  Carlo  fluctuations.  This  behavior  is  ex¬ 
pected,  based  on  the  preceding  discussion.  However, 
we  also  note  that  all  models  predict  nearly  the  same 
values  for  Ed  and  Eou,  which  is  counter  to  intuition 
based  on  oceanographic  experience.  This  result  is 
easily  explained  if  we  recall  that  the  Rayleigh  phase 
function  is  nearly  isotropic  (independent  of  the  scat¬ 
tering  angle)  and  that  the  medium  is  highly  scattering. 
Because  of  the  intense  scattering,  the  incident  colli¬ 
mated  radiance  distribution  approaches  its  asymptotic 
form  very  quickly  with  depth.  Preisendorfer49  shows 
that  for  an  isotropic  phase  function  the  asymptotic 
radiance  distribution,  L«,  has  an  elliptical  shape: 


LJd)  = 


Lp 

1  +  cos  0 


(10) 


Here  L0  depends  only  on  the  inherent  optical  proper¬ 
ties  and  k„  is  the  eccentricity  of  the  ellipse;  is 
numerically  equal  to  the  nondimensional  asymptotic 
diffuse  attenuation  coefficient.  The  analytic  forms 
of  Loo  for  a  Rayleigh  phase  function  and  a  Rayleigh 
phase  matrix  are  also  known.50  For  w0  =  0.9  the 
Rayleigh  L„  is  very  close  to  elliptical,  and  so  we  can 
use  the  simpler  form  of  Eq.  (10)  for  the  following 
argument.  The  Ed  and  Eou  corresponding  to  L„  of 
Eq.  (10)  are 


ZTTlsn  r 

Ed - —[k.  +  lnd-ku,)], 

’*00 

2irL0 

Eou  =  -~ln(l  +  *.).  (11) 

Now  the  value  of  for  the  problem  at  hand  turns  out 
to  be  ku,  =  0.52  (see  Table  7).  This  value  is  coinciden¬ 
tally  very  near  to  the  value  kx  =  0.531,  which  makes 
Ed  =  Eou  in  Eq.  (11),  thus  explaining  the  numerical 
results  seen  in  Fig.  2(a).  This  peculiar  behavior  of 
Ed  and  Eou  depends  on  both  the  phase  function  and 
the  scattering-to-attenuation  ratio.  Such  behavior 
is  not  seen  in  the  output  for  the  other  problems,  nor 
would  it  ever  be  encountered  in  a  natural  water  body. 

Note  also  that  both  Ed  and  Eou  are  greater  just 
below  the  water  surface  than  just  above  it,  which  may 
also  seem  counterintuitive.  However,  this  is  just  the 
phenomenon  of  optical  energy  trapping  in  highly 
scattering  waters,  as  discussed  by  Stavn  et  a/.51  and 
by  Plass  et  al .52  In  the  present  case  of  a  solar  angle 
of  60°,  more  than  93%  of  the  incident  solar  irradiance 
is  transmitted  through  the  level  surface  into  the 
water.  Approximately  one  half  of  the  highly  diffuse 


20  December  1993  /  Vol.  32,  No.  36  /  APPLIED  OPTICS  7495 


upwelling  irradiance  just  below  the  surface  is  re¬ 
flected  back  down  by  the  surface.  The  total  Ed  just 
below  the  surface  is  the  sum  of  the  transmitted  solar 
contribution  and  the  reflected  upwelling  contribu¬ 
tion;  this  sum  is  greater  than  Ed(air).  Likewise, 
£ou(air)  consists  of  the  (relatively  weak)  specularly 
reflected  solar  beam  plus  diffuse  light  transmitted 
upward  through  the  water  surface;  this  sum  is  less 
than  Eou  just  below  the  surface. 

Figure  2(b)  shows  the  output  for  the  Rayleigh 
phase  function  and  a  highly  absorbing  medium  with 
wo  =  0.2.  Now  Eou  is  an  order  of  magnitude  less  than 
Ed.  There  is  a  spread  of  almost  a  factor  of  3  in  the 
Monte  Carlo  estimates  of  Eou  at  t  =  10,  and  three  of 
the  Monte  Carlo  models  had  too  few  photons  left  at 
t  =  10  to  provide  an  estimate  of  Lu  at  that  depth. 
This  behavior  is  expected  for  this  highly  absorbing 
case. 

Table  5  displays  the  average  (over  all  models) 
values  of  Ed<  Eou,  and  Lu  at  selected  depths  for  this 
and  the  remaining  problems.  These  data  are  pro¬ 
vided  for  readers  who  wish  to  compare  their  own 
models  with  ours.  Such  comparisons  should  be  espe- 
cialy  worthwhile  for  simple  parameterized  models 
that  attempt  to  compute  irradiances  without  solving 
the  complete  radiative  transfer  equation.  Table  5 
also  displays  the  ratio  of  the  sample  standard  devia¬ 
tion  s  to  the  sample  mean  x, 


s 

x 


N 


N 


—  X  (*i  -  xf 

“  A  i=l 


1/2 


1  N 


where  x,  is  the  result  predicted  by  the  ith  model  for 
the  quantity  of  interest  and  N  is  the  number  of  model 
predictions  (N  =  7  for  most  quantities).  The  ratio  s/x 
is  a  quantitative  measure  of  how  close  together  the 
models’  predictions  are  for  a  given  quantity.  Inspec¬ 
tion  of  this  ratio  for  problem  1  shows  that  the  model 
predictions  are  usually  closer  together  for  the  highly 
scattering  case  (w0  =  0.9)  than  for  the  highly  absorb¬ 
ing  case  (wo  =  0.2),  closer  together  at  shallow  depths, 
and  closest  together  for  Ed.  The  greatest  spread  in 
values  is  for  Lu  at  large  depths,  because  of  the  small 
number  of  photons  available  for  its  estimation  by  the 
Monte  Carlo  models. 

Problem  2.  Figure  3  shows  the  models’  output  for 
problem  2.  Figure  3(a)  is  for  the  highly  scattering 
case  of  to0  =  0.9.  Each  of  the  seven  models  provides 
essentially  the  same  values  for  Ed  and  for  Eou  to  10 
optical  depths  (and  deeper);  some  Monte  Carlo  fluctua¬ 
tion  is  apparent  in  the  Lu  values.  Figure  3(b)  shows 
the  same  computations  for  the  highly  absorbing  case 
of  to0  =  0.2.  Once  again,  all  models  give  nearly  the 
same  values  for  Ed  and  for  Eou  to  10  optical  depths. 
Now,  however,  considerable  Monte  Carlo  fluctuation 
in  the  Lu  values  is  seen  at  even  shallow  depths;  only 
models  II,  DO,  and  MC3  were  able  to  compute  Lu 
below  t  =  10. 


We  emphasize  that  the  large  fluctuations  seen  in 
some  of  the  estimates  in  Fig.  3(b)  are  simply  the 
result  of  tracing  an  insufficient  number  of  photons  in 
the  simulations,  and  not  of  any  inadequacies  in  the 
models  themselves.  Tracing  additional  photons,  at  a 
proportional  increase  in  computational  expense,  can 
reduce  these  fluctuations  to  any  desired  level.  The 
particular  values  seen  in  Fig.  3  are  each  the  result  of 
one  simulation.  Running  the  Monte  Carlo  models 
with  different  seeds  for  their  random  number  genera¬ 
tors  would  generate  a  noticeably  different  set  of 
curves  for  those  instances  where  large  fluctuations 
are  seen  in  Fig.  3.  It  should  be  noted  that  there  are 
certain  sampling  schemes  that  can  improve  the  statis¬ 
tics  at  greater  depths.  However,  this  improvement 
is  usually  at  the  expense  of  larger  errors  in  the 
radiometric  quantities  at  smaller  depths. 

The  euphotic  zone  is  the  region  of  a  water  body 


Table  5.  Average  Values  of  and  l*  at  Selected  Depths  for 
Problems  1-6° 


Optical 

Depth 

Average  Value 

Corresponding  s  /x 

Ed  Eou 

Lu 

Ed 

Eou 

Lu 

Problem  1,  wo  =  0.9  (N  = 

7) 

1 

3.66-1  3.72-1 

4.85-2 

0.002 

0.005 

0.015 

5 

4.33-2  4.35-2 

5.59-3 

0.003 

0.007 

0.052 

10 

3.16-3  3.20-3 

4.37-4 

0.015 

0.038 

0.091 

Problem  1,  wo  =  0.2  ( N  = 

7) 

1 

1.41-1  1.34-2 

1.72-3 

0.001 

0.003 

0.044 

5 

1.07-3  1.00-4 

1.37-5 

0.005 

0.039 

0.288 

10 

2.93-6  3.00-7 

3.39-8  (N  =  4) 

0.102 

0.308 

0.197 

Problem  2,  u>o  =  0.9  (JV  = 

7) 

1 

4.13-1  9.31-2 

6.99-3 

0.001 

0.021 

0.063 

5 

1.87-1  4.63-2 

3.26-3 

0.005 

0.017 

0.055 

10 

6.85-2  1.65-2 

1.21-3 

0.010 

0.014 

0.109 

Problem  2,  wo  =  0.2  ( N  = 

7) 

1 

1.62-1  9.66-4 

5.47-5 

0.000 

0.023 

0.060 

5 

2.27-3  1.37-5 

6.24-7  (JV  =  6) 

0.002 

0.063 

0.355 

10 

1.30-5  7.28-8 

4.02-9  (N  =  5) 

0.047 

0.187 

0.248 

*  Problem  3  (N  =  6) 

5  m 

2.30-1  4.34-2 

3.13-3 

0.006 

0.025 

0.054 

25  m 

1.62-3  2.86-4 

2.12-5 

0.028 

0.038 

0.061 

60  m 

5.23-5  5.13-6 

3.57-7 

0.071 

0.036 

0.434 

Problem  4  (IV  =  6'f’ 

1 

3.23-1  7.13-2 

5.63-3 

0.076 

0.091 

0.111 

5 

1.49-1  3.57-2 

2.77-3 

0.072 

0.076 

0.141 

10 

5.56-2  1.31-2 

9.60-4 

0.070 

0.073 

0.107 

Problem  5  (JV  =  4) 

1 

1.14-1  3.55-2 

2.09-3 

0.012 

0.020 

0.031 

5 

4.33-2  1.22-2 

7.63-4 

0.009 

0.028 

0.036 

10 

1.48-2  3.65-3 

2.49-4 

0.007 

0.020 

0.025 

Problem  6  (N  =  3) 

1 

1.62-1  9.81-4 

6.84-5 

0.000 

0.010 

0.020 

5 

2.28-3  2.28-3 

3.60-4 

0.003 

0.002 

0.010 

“N  is  the  number  of  models  included  in  the  averages.  The  ratio 
of  the  sample  standard  deviation  to  the  sample  mean,  s/x,  is  also 
displayed  for  each  average  value.  The  average  values  are  relative 
to  an  incident  solar  irradiance  of  EL  =  1.0  W  m"2  nm"1  upon  the 
water  surface,  except  for  problem  4,  for  which  Ex  is  applied  at  the 
top  of  the  atmosphere.  The  notation  3.66- 1,  for  example,  means 
3.66  x  10-'. 

bs/x  values  determined  by  systematic  offset;  see  discussion  in  the 
text. 
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(a) 


(b) 

Fig.  3.  Model  predictions  for  problem  2,  the  base  case:  (aj  u>o  = 
0.9  and  (b)  o>o  =  0.2. 


where  there  is  sufficient  light  for  photosynthesis  to 
take  place.  In  normal  daylight  conditions,  it  extends 
from  the  surface  to  a  depth  where  the  irradiance  is 
roughly  1%  of  its  surface  value.  We  see  in  Fig.  3(b) 
that  Ed  and  Eou  have  decreased  by  2  orders  of 
magnitude  at  approximately  4  optical  depths.  Each 
of  the  models  produces  nearly  identical  irradiances  to 
depths  greater  than  t  =  4,  so  that  each  of  the  models 
is  perfectly  adequate  for  the  purposes  of  biological 
oceanography.  Likewise,  the  models  produce  very 
nearly  the  same  water-leaving  radiances,  Lu( air),  as 
would  be  of  interest  in  remote-sensing  studies. 

Problem  3.  Figure  4  shows  the  models’  output  for 
problem  3,  the  stratified  water  case.  The  1%  irradi¬ 
ance  level  is  now  at  approximately  z  =  21  m.  Once 
again,  the  models  provide  nearly  identical  output  to 
depths  far  below  the  euphotic  zone. 

Problem  4.  Figure  5  shows  Ed  values  near  the 
water  surface  for  the  simulation  of  problem  4,  the 
case  with  an  atmosphere.  The  different  ways  in 
which  the  models  simulate  the  atmosphere  lead  to  an 
18%  spread  in  the  values  of  Ed  just  above  the  water 


Fig.  4.  Model  predictions  for  problem  3,  the  stratified-water  case. 


surface.  This  difference  in  Ed(air)  values  is  then 
carried  throughout  the  underwater  computations. 
The  s/x  ratio  displayed  in  Table  5  is  uniformly  large 
for  this  problem  because  of  the  systematic  offset  of 
the  different  models’  predictions.  Note  that  appar¬ 
ent  optical  properties,  such  as  reflectances  and  diffuse 
attenuation  functions,  are  not  affected  by  this  offset, 
because  the  apparent  properties  are  defined  as  ratios 
of  radiometric  quantities.  For  example,  the  s /x.  ratio 
for  the  Kd  values  computed  from  the  plotted  Ed  values 
at  depths  z  =  0  and  1  m  is  0.009,  which  is  much 
smaller  than  the  s/x  =  0.076  value  tabulated  for  Ed 
atT  =  1. 

Problem  5.  Four  of  the  models  (II,  MCI,  MC3,  and 
MC4)  are  capable  of  simulating  a  windblown  air- 
water  surface  as  defined  in  problem  5.  Figure  6 
shows  output  from  these  models  for  a  solar  zenith 
angle  of  0^  =  80°.  The  models  are  nearly  identical 
in  their  output,  even  in  this  case  of  nearly  horizontal 
incidence,  for  which  any  differences  in  the  models 
should  be  most  noticeable.  Note  that  Z?ou(air)  is 


Fig.  5.  Ed  near  the  surface  for  problem  4,  the  base  case  plus  an 
atmosphere. 


20  December  1993  /  Vot.  32,  No.  36  /  APPLIED  OPTICS  7497 


Fig.  6.  Model  predictions  near  the  surface  for  problem  5,  the 
capillary-wave  case.  The  wind  speed  is  U  =  7.23  m  s-1,  and  the 
zenith  angle  of  the  sun  is  03U„  =  80°. 


greater  than  Ed(air).  This  is  because  Eou( air)  con¬ 
tains  a  large  contribution  by  the  specularly  reflected 
solar  beam:  simulations  by  Preisendorfer  and  Mob¬ 
ley4  show  that  the  reflectance  of  a  capillary-wave 
surface  is  greater  than  0.22  for  a  wind  speed  of  7.23  m 
s’1  and  0sun  =  80°.  The  solar  beam  contribution  to 
Ea  is  weighted  by  a  cos  08un  factor,  which  is  small  for 
0SU„  =  80°. 

Problem  6.  Models  II,  DO,  and  MC3  can  simulate 
a  finite-depth  bottom.  Figure  7  shows  the  output 
from  both  models  for  the  case  of  o>0  =  0.2;  the  models 
are  clearly  in  excellent  agreement.  It  is  easy  to  show 
that  Eou  =  Ed  for  a  Lambertian  surface  of  reflectance 
0.5,  and  all  three  models  show  this  expected  result  at 
depth  t  =  5. 

Problem  7.  Four  of  the  models  (MCI,  MC2,  MC3, 
and  MC5)  can  simulate  Raman  scattering.  Table  6 
compares  the  inelastically  scattered  contributions  to 
the  downwelling  and  upwelling  plane  irradiances,  Ed 
and  Eu,  respectively,  for  the  simulation  defined  in 


in  air 
0.0 

h  1.0 

x 

CL 

« 

^  2.0 

o 

u 

g-  3.0 


4.0 
5.0 

10-5  10-4  10~3  10~2  10'1  10° 
spectral  irradlance  E  (W  m"'  nm'1) 
spectral  radiance  Lu  (W  m"  sr"  nm") 


Fig.  7.  Model  predictions  for  problem  6,  the  finite-depth  case. 
The  bottom  reflectance  is  0.5. 


Table  6.  Raman  Scattering  Contribution*  to  £„  and  £„  at  g  =  4S6  nm 
From  an  Excitation  Wavelength  of  c  „  =  417  nm* 


Depth 

(m) 

Model 

MCI 

MC2 

MC3 

MC5 

0 

0.01875 

Ed  values 
0.01874 

0.01739 

0.01873 

50 

0.02489 

0.02488 

0.02470 

0.02490 

100 

0.01136 

0.01136 

0.01123 

0.01138 

0 

0.03532 

Eu  values 
0.03512 

0.03478 

0.03523 

50 

0.01034 

0.01042 

0.01027 

0.01039 

100 

0.00287 

0.00296 

0.00292 

0.00296 

“Parameter  values  are  given  in  the  specification  of  problem 
7.  Values  in  the  body  of  the  table  have  units  of  W  m~2  nm1  for  an 
incident  irradiance  ofJSi  =  1.0  W  m'2  nm-1  at 


problem  7 .  The  models  are  clearly  in  excellent  agree¬ 
ment,  even  though  their  respective  formulations  of 
inelastic  scatter  are  somewhat  different. 

Computation  of  radiance  distributions.  Five  of 
the  models  (II,  DO,  MC2,  MC3,  and  MC4)  compute 
the  full  radiance  distribution,  rather  than  just  tally¬ 
ing  photons  as  necessary  to  compute  the  irradiances 
and  Lu.  Figure  8  illustrates  the  consistency  with 
which  the  various  models  compute  the  radiance  distri¬ 
bution.  The  figure  shows  L(t,  0,  4>)  in  the  plane  of 
the  sun  at  depths  of  t  =  0,  5,  and  20  for  probl' m  2, 
wo  =  0.9.  Direction  (0„,  (j>„)  gives  the  viewing  direc¬ 
tion,  i.e.,  the  direction  an  instrument  points  to  detect 
photons  traveling  in  the  (0  =  180°  -  0„,  <}>  =  180°  +  <J>„) 
direction.  Thus  0„  =  180°  corresponds  to  looking 
straight  up  and  seeing  photons  heading  straight 
down;  the  nadir  radiance,  Lu,  of  Fig.  3(a)  is  the  value 
plotted  at  0„  =  0°.  The  sun  is  in  the  <)>„  =  0° 
half-plane. 

The  curves  of  Fig.  8  are  explained  as  follows.  We 
begin  at  t  =  0  (in  the  water  just  below  the  surface) 


viewing  direction  (ev,py)  (deg) 

Fig.  8.  Radiance  distribution  in  the  plane  of  the  sun  for  problem 
2,  wo  =  0.9.  Angles  (0„  <M  are  viewing  directions:  0„  =  180 0  -  0 
and  <t>„  =  180°  +  <t>,  where  (0,  <j>)  are  the  directions  of  photon 
travel.  The  solid  curves  are  L( t,  0„,  <j>„)  at  selected  depths  within 
the  water  for  models  II  and  DO;  models  MC2-MC4  are  shown  by 
the  dashed  curves.  The  dotted  curve  is  the  asymptotic  dis 
tion  L*(0„)  normalized  to  the  largest  value  of  L  at  t  =  20. 
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with  our  backs  to  the  sun  (looking  in  the  «J>t.  =  180° 
direction).  Looking  straight  down  we  see  the  nadir 
radiance  at  (8„,  4>L, )  =  (0°,  180°).  Looking  up  toward 
the  horizontal  (0U  =  90°),  the  radiance  increases 
slightly  because  of  total  internal  reflection  of  radiance 
that  has  been  scattered  into  nearly  horizontal  direc¬ 
tions.  The  radiance  then  decreases  quickly  as  our 
viewing  angle  passes  beyond  the  critical  angle  for 
total  internal  reflection.  In  the  region  around  0„  = 
180°  we  are  looking  upward  and  seeing  the  upwelling 
radiance  that  is  reflected  downward  by  the  level  water 
surface.  Note  for  example  (using  the  digital  output 
from  Model  II)  that  L(t  =  0,  0„  =  180°)/L(t  =  0, 
0„  =  0°)  =  1.737  x  10-4/8.236  x  10~3  =  0.021,  which 
is  just  the  Fresnel  reflectance  of  the  surface  for 
perpendicular  incidence.  Recall  that  in  problem  2 
the  sky  is  black,  so  there  is  no  sky  radiance  transmit¬ 
ted  through  the  surface.  In  problem  4  (not  shown), 
transmitted  sky  radiance  fills  in  the  large  dip  in  the 
radiance  near  8„  =  180°.  As  our  view  passes  the 
zenith  we  are  now  facing  the  sun.  The  large  spike  in 
the  radiance  near  (0„,  4>J  =  (140°,  0°)  is  the  refracted 
solar  beam.  The  noticeable  0„  offset  in  the  position 
of  the  plotted  peak  radiance  occurs  because  different 
models  choose  their  quad  boundaries  differently. 
The  radiance  values  are  plotted  at  the  0„  values  of  the 
quad  centers,  which  range  from  135°  to  139.7°  for  the 
quad  containing  the  refracted  solar  beam;  plotted 
points  are  connected  by  straight  lines.  Looking  be¬ 
yond  the  sun,  we  see  a  large  horizontal  radiance, 
which  decreases  as  we  look  downward. 

Model  DO  shows  a  more  pronounced  spike  in  the 
radiance  near  the  solar  direction,  and  more  pro¬ 
nounced  changes  near  the  critical  angle  than  do  the 
other  models.  This  is  because  model  DO  computes 
radiances  in  specific  directions,  rather  than  quad- 
averaged  radiances.  The  angular  quadrature  points 
in  model  DO  are  clustered  near  the  critical  angle  and 
near  the  horizon,  to  get  increased  resolution  in 
regions  where  the  radiance  varies  rapidly  with  polar 
angle. 

By  depth  r  =  5,  scattering  has  smeared  out  the 
solar  beam  and  increased  the  downwelling  radiance 
seen  when  looking  upward  near  the  zenith.  The 
radiance  distribution  at  t  =  20  is  very  similar  in  shape 
to  the  asymptotic  distribution,  L®(0„).  The  asymp¬ 
totic  distribution  as  computed  by  model  II  and  normal¬ 
ized  to  the  largest  value  of  L(r  =  20)  is  shown  as  a 
dotted  line  in  Fig.  8.  Note  that  only  a  small  amount 
of  Monte  Carlo  fluctuation  is  seen  even  at  t  =  20,  for 
this  highly  scattering  case. 

Radiance  distributions  computed  by  the  various 
models  are  in  equally  close  agreement  for  the  other 
canonical  problems  (except  for  Monte  Carlo  fluctua¬ 
tions  in  the  small  o>0  cases)  and  will  not  be  discussed. 

Computation  of  asymptotic  radiances.  The  asymp¬ 
totic  radiance  regime  (also  called  the  diffusion  regime) 
is  the  region  far  enough  from  the  boundaries  of  a 
homogeneous  medium  that  the  radiance  is  indepen¬ 
dent  of  the  incident  direction  of  the  source  and  of 
boundary  effects.  Radiance  in  the  asymptotic  re¬ 


viewing  direction  0,  (deg) 


Fig.  9.  Asymptotic  radiance  distributions  LJ8r)  for  problems  1 
and  2,  as  computed  by  various  models  (solid  curves).  The  dotted 
curves  give  the  exact  analytic  solution46  for  the  Rayleigh  phase 
function  of  problem  1. 

gime  is  independent  of  the  azimuthal  angle  and  it 
decreases  exponentially  as  exp(-fc*T).  Here  kx  = 
Kx/c,  where  if*  is  the  dimensional  (in  inverse  meters) 
asymptotic  diffuse  attenuation  coefficient  and  c  is  the 
beam  attenuation  coefficient.  The  shape  LJB)  of  the 
asymptotic  radiance  distribution  is  determined  only 
by  the  inherent  optical  properties  of  the  water;  it  is 
independent  of  depth. 

Model  II  computes  L„(0)  and  the  associated  value  of 
the  nondimensional  asymptotic  diffuse  attenuation 
coefficient  kx  by  the  solution  of  a  matrix  eigenvalue 
equation.6  The  smallest  eigenvalue  of  the  matrix  is 
kx,  and  the  associated  eigenvector  gives  Lx.  Model 
DO  obtains  the  asymptotic  solution  in  a  similar 
fashion.  Models  MC2  and  MC3  obtain  Lx  and  kx  by 
solution  of  the  equivalent  integral  equation49  53 

J»2ir  /»1 

Lx(p/)P(»|i)dp'd<}>'.  (12) 

o  J- 1 

The  exact  analytical  solution  to  Eq.  ( 12)  for  the  case  of 
scattering  according  to  a  Rayleigh  phase  function,  as 
well  as  for  a  Rayleigh  phase  matrix,  was  found  by 
Kattawar  and  Plass.50  Numerical  solutions  for  phase 
functions  that  are  highly  peaked  in  the  forward 
direction  have  been  given  in  Kattawar  and  Plass50 
and  in  Prieur  and  Morel.54 

Figure  9  shows  the  computed  Lx(9j,  normalized  to 
one  at  0„  =  180°,  for  problems  1  and  2.  Table  7 


Table  7.  Computed  Values  of  k 


Problem 

Model 

ii 

DO 

MC1“ 

MC2 

MC3 

1,  t»o  -  0.9 

0.5248 

0.5232 

0.52 

0.5232 

0.5235 

1,  (i>0  =  0.2 

1.0006 

0.9994 

— 

0.9996 

0.9952 

2,  ii>o  —  0.9 

0.1920 

0.2068 

0.189 

0.1835 

0.1879 

2,  ioq  ~  0.2 

0.8737 

0.8794 

— 

0.8590 

0.8619 

“Values  determined  by  visual  inspection  of  plotted  output. 
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Table  8.  Comparison  of  Percent  Accuracies  for  Computing  and 
Measuring  Radiometric  Variables 


Variable 

2a  Spread 
of  Model 
Values 

Current 

Measurement 

Capability 

Target 
Accuracy 
for  SeaWiFS“ 

Ed 

1 

3-5 

2 

Eou 

5 

3-5 

— 

Lu 

12 

3-5 

3 

“From  Mueller  and  Austin.55 


shows  the  corresponding  A„  values.  The  numerical 
results  are  in  excellent  agreement  for  problem  2  and 
for  the  wo  =  0.9  case  of  problem  1,  which  also  agrees 
with  its  exact  analytic  solution.  However,  the  nu¬ 
merical  results  differ  considerably  for  the  wo  =  0.2 
case  of  problem  1,  and  each  is  considerably  off  from 
the  analytic  solution.  The  reason  for  this  inaccuracy 
in  the  computed  L„  is  as  follows.  For  problem  1, 
o)0  =  0.2,  the  analytic  A„  value  is  A.  =  0.99937. 
However,  Eq.  (12)  becomes  singular  as  p.  -»  1  when 
A*  =  1.  For  the  nearly  singular  case  at  hand,  both 
model  II’s  eigenmatrix  routine  and  the  integral  equa¬ 
tion  routines  are  having  a  difficult  time  determining 
accurate  values  for  A„  and  L„.  This  is  most  notice¬ 
able  in  the  A„  =  1.0006  value  determined  by  model  II; 
the  theoretical  upper  limit  for  A»  is  exactly  one. 
Even  slight  errors  in  A„  cause  large  differences  in  L* 
when  A.  is  near  one.  Kattawar  was  able  to  obtain  a 
satisfactory  numerical  solution  of  Eq.  (12)  for  this 
case  only  after  resorting  to  quadruple-precision  arith¬ 
metic.  The  A„  *  0.87  value  seen  in  problem  2,  o)0  = 
0.2,  is  far  enough  from  one  that  no  numerical  difficul¬ 
ties  arise.  Note  that  the  computation  of  A.  and  L„  is 
a  separate  problem  from  the  computation  of  the 
radiances  and  irradiances  as  discussed  above.  The 
inaccuracies  in  A„  and  just  discussed  in  no  way 
imply  inaccuracies  in  the  solution  of  Eq.  ( 1). 

5.  Conclusions 

Problems  1-3  of  Section  3  cover  the  extreme  range  of 
oceanic  inherent  optical  properties:  u>0  from  0.2  to 
0.9,  phase  functions  for  pure  Rayleigh  and  pure 
particulate  scattering,  and  strong  vertical  stratifica¬ 
tion.  In  computations  of  Ed  and  Eou,  the  numerical 
models  of  Section  2  usually  gave  results  within  a  few 
percent  of  each  other  throughout  the  euphotic  zone. 
The  spread  in  Lu  values  was  as  large  as  12%  in  highly 
scattering  waters  and  much  larger  in  highly  absorb¬ 
ing  waters  at  the  bottom  of  the  euphotic  zone. 

The  statistical  fluctuations  of  the  Monte  Carlo 
results  from  the  true  values  of  the  predicted  quanti¬ 
ties  are  normally  distributed.  We  therefore  expect 
that  more  than  95%  of  the  Monte  Carlo  simulations 
will  be  within  2  standard  deviations  (2oj  of  the  correct 
value.  The  data  of  Table  5  give  us  a  feeling  for  the 
size  of  this  2cr  spread  of  values.  Table  8  shows  the  2ct 
spread  (expressed  as  a  percentage  of  the  mean)  for  Ed , 
E^  and  Lu  in  near-surface  waters  (based  on  t  =  1  for 
problems  1  and  2  and  based  on  z  =  5  m  for  problem  3). 
Column  2  of  Table  8  shows  typical  errors  in  these 


radiometric  quantities  when  measured  by  commer¬ 
cial  instruments  now  in  wide  use.  The  third  column 
of  the  table  shows  the  accuracy  desired  in  measure¬ 
ments  to  be  used  for  ground-truth  validation  of  the 
SeaWiFS  ocean  color  satellite55  (to  be  launched  in 
1994).  Obtaining  such  accuracies  in  Ed  and  Lu  mea¬ 
surements  requires  very  careful  instrument  calibra¬ 
tion. 

We  see  from  Table  8  that  the  present  numerical 
models  easily  compute  Ed  with  greater  accuracy  than 
can  be  obtained  with  current  instruments.  Numeri¬ 
cal  estimates  of  Eou  have  approximately  the  same 
accuracy  as  measured  values.  The  computed  values 
of  Lu  are  less  accurate  than  can  be  measured  or  than 
are  needed  for  remote-sensing  studies  requiring  abso¬ 
lute  radiometric  values  of  Lu .  Thus  the  Monte  Carlo 
models  should  trace  more  photon  histories,  if  very 
accurate  Lu  values  are  required.  The  standard  devia¬ 
tion  of  the  Monte  Carlo  fluctuations  is  proportional  to 
n'1/2,  where  n  is  the  number  of  photons  traced. 
Therefore  the  2cr  spread  seen  in  Table  8  can  be  cut  in 
half  by  tracing  four  times  as  many  photons,  which  is 
computationally  practicable.  Another  possibility  is 
to  use  the  backward  Monte  Carlo  method,  as  de¬ 
scribed  in  Gordon.56 

Monte  Carlo  calculations  made  using  statistical 
estimation  techniques  can  also  yield  continuum  radi¬ 
ances,  rather  than  quad-averaged  values.  Thus  if 
one  is  interested  in  results  for  a  few  detectors  located 
at  precise  angles,  this  technique  can  give  highly 
accurate  radiance  values  with  only  a  very  few  photons 
being  traced.57"59 

Values  predicted  by  the  Monte  Carlo  models  gener¬ 
ally  fall  on  both  sides  of  the  values  predicted  by 
models  II  and  DO,  which  do  not  have  statistical 
fluctuations.  Thus  models  II  and  DO  have  an  advan¬ 
tage  in  the  computation  of  upwelling  quantities  or  in 
computations  at  great  depths,  which  require  tracing 
very  large  numbers  of  photons  in  the  Monte  Carlo 
codes. 

The  systematic  differences  in  the  atmospheric  mod¬ 
els  used  to  simulate  problem  4  lead  to  a  2o  spread  of 
the  order  of  20%  in  the  computed  radiometric  quanti¬ 
ties.  Thus  to  compute  acceptably  accurate  absolute 
radiometric  values,  more  careful  attention  must  be 
paid  to  how  the  incident  radiance  upon  the  water 
surface  is  obtained.  However,  as  noted  before,  sys¬ 
tematic  offsets  in  the  absolute  radiometric  variables 
do  not  affect  the  values  of  apparent  optical  properties 
obtained  from  the  radiometric  variables.  The  pres¬ 
ent  simple  atmospheric  models  therefore  all  appear  to 
be  satisfactory  for  the  computation  of  apparent  opti¬ 
cal  properties. 

Based  on  the  problem  solutions  presented  above, 
and  on  such  comparisons  between  models  and  oceano¬ 
graphic  measurements  as  have  been  made  (not  dis¬ 
cussed  here),  we  conclude  that  each  of  the  numerical 
models  discussed  here  incorporates  correct  math¬ 
ematical  representations  of  the  relevant  radiative 
processes  (absorption  and  elastic  and  inelastic  scatter¬ 
ing)  and  of  the  effects  of  the  air-water  boundary. 


7500  APPLIED  OPTICS  /  Vol.  32,  No.  36  /  20  December  1993 


Moreover,  the  models  provide  accurate  numerical 
solutions  of  the  associated  equations.  Each  of  these 
models  is  adequate  for  most  of  the  needs  of  optical 
oceanography  and  limnology. 

Appendix  A:  inelastic  Source  Function  for  Model  MCI 

As  noted  in  Section  2,  model  MCI  incorporates 
Raman  scattering  (and  other  inelastic  processes,  such 
as  fluorescence)  in  an  azimuthally  averaged  form 
suitable  for  the  computation  of  inelastic-scattering 
effects  on  irradiances.  The  corresponding  math¬ 
ematical  form  of  the  source  function,  which  is  used  in 
the  4>-averaged  version  of  Eq.  (1),  is  developed  as 
follows.  This  formulation  is  based  on  expanding 
both  the  Raman  scattering  function  and  the  azimu- 
multally  averaged  excitation  radiance  in  a  series  of 
Legendre  polynomials.60 

The  source  function  for  inelastic  processes  is  given 

by 


sjz,  e,  x) 

i  *>  r 

»  2  I  binl)(2<  Kx  -*  X)P,( COS  0)£/(.Z,  Xex)dXex, 
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gendre  polynomials  and  the  quantity  y  s  (1  -  p) 
(1  +  3p)  gives 
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Comparing  Eqs.  ( Al )  and  (A3 )  reveals  that 

^Ram  (*>  Xex  *  X)  Xex  *  X), 

1  ~  P 

1  +  2p 


bKam(Z,  K 


X), 


and  that  all  of  the  other  are  zero.  Finally,  the 
4>-averaged  source  function  for  Raman  scattering  is 
given  by 


Snam(z 


bRam(*j  Xe: 


X)EoU,  Xex) 
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1/  1-p  \E2(z,Kx) 
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In  general,  at  the  emission  wavelength  X,  the 
source  function  resulting  from  a  narrow  band  of 
excitation  wavelengths  AXex  is 


Sjz,  0,  X)  =  ^  bjz,  X, 


X)AxcxFo(2,  Xex) 
"  bj\z,  Kx  X)jE/(z,  Xex) 


X  2 

/  =  0  bJz,  Kx  *  K)Eq(z,  Kx) 


Pl( COS  0). 
(A4) 


In  these  equations,  Pi  is  the  Legendre  polynomial  of 
order  l,  N  =  0  for  isotropically  emitted  fluorescence, 
and  IV  =  2  for  Raman  scattering.  The  total  inelastic¬ 
scattering  coefficient  6,n(0)  is 


bj0)(z,  Kx  -*  X)  3  bjz,  Xex  K 

•  J  flin(z;  0’,  4)’ 


0, 4>;  Kx  —  x)dfl’ 


=  2tt  flin(z,  ||/,  Xex  -*  X)sin  t|idi|i. 

Jo 


Ei  for  /  =  0  is  just  the  scalar  irradiance  at  Xex,  whereas 
Ei  for  l  =  1  is  the  net  irradiance  Ed  -  Eu  at  Xex.  The 
inelastic  component  of  the  irradiance  at  X  depends 
only  on  the  irradiances  at  the  excitation  wavelength(s) 
and  on  the  bj!)  coefficients  for  the  particular  process. 

For  Raman  scattering,  flin  s  0^,  and  the  angular 
distribution  of  0^  is  given  by 

flRam(*,  X«  X)  =  flRam(*l»)&Ram(*>  Kx  “ *  X),  (A2) 

where  flitamW  is  given  by  Eq.  (9).  Substituting  Eq. 
(9)  into  Eq.  (A2)  and  rewriting  in  terms  of  the  Le- 


To  simulate  the  irradiances  at  X,  the  basic  Monte 
Carlo  code  for  elastic  scattering  only  is  run  at  Xex  to 
determine  Ei(z,  Xex).  Then  Eq.  (A4)  is  used  to  inject 
inelastically  scattered  photons  into  the  medium  with 
the  proper  distribution  in  z  and  0.  One  way  to  do 
this  is  to  choose  z  from  the  probability  density  p(z) 
given  by 


p(z)  = 


Eq(z,  K 


f  E0(z, 

Jo 


Xex)d Z 


Thus,  given  a  random  number  p,  from  the  sequence 
■  •  ■  Pj,  Py+i,  p>+2,  •  •  - ,  z  is  found  from 


>-r 

*0 


p(z')dz'. 


Given  z,  we  then  choose  0  from  the  conditional 
density  p(djz)  given  by 


p(0|z) 


4l T 


1 

/=0 


bin i!'(z,  Kx  ->  X)E;(Z,  Xex) 
b,n (z,  Xex  *  X)E„(z,  Xex) 


Pi(cos  0), 
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so  that 

Py  +  l  =  [  p(O'lz)d0’. 

Jo 

Finally,  the  rest  of  the  source  function  must  be 
incorporated  into  a  photon  weight 

Ax 

W  =  bin(z,  Xex  -*  X)AXex  E0{z,  X  Jdz, 

Jo 

so  that  Sin(z,  0,  X)  =  Wp(z)p(0 1 z)  as  required.  Once  a 
photon  is  emitted  at  X,  it  is  followed  in  a  manner 
similar  to  that  which  would  be  used  in  the  absence  of 
inelastic  scattering,  with  the  exception  that  inelastic 
scattering  from  X  to  longer  wavelengths  is  included  by 
increasing  the  absorption  coefficient  a(z,  X)  by  the  ap¬ 
propriate  inelastic-scattering  coefficient  bm(z,  X  -*  X') 
with  X  <  X'  (recall,  however,  the  discussion  of  this 
point  in  the  definition  of  problem  7).  In  the  code, 
E0(z,  Xex)  is  normalized  to  unit  irradiance  at  X«x  enter¬ 
ing  the  top  of  the  atmosphere  normal  to  the  solar 
beam,  so  the  computed  irradiances  at  X  (as  seen  in 
Table  6)  are  for  unit  irradiance  at  Xex  entering  the  top 
of  the  atmosphere.  This  simulation  technique  was 
in  part  developed  in  this  manner  so  that  it  could  be 
used  with  experimental  measurements  of  Ei  to  pre¬ 
dict  the  inelastically  scattered  irradiances  for  a  given 
process,  e.g.,  Raman  scattering.  Such  measure¬ 
ments  can  be  obtained  using  instrumentation  devel¬ 
oped  by  Voss.61-62 
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